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Foreword 

One  of  the  most  famous  American  physicists  of  the 
twentieth  century,  Richard  Feynman,  in  1982  was  the 
first  to  propose  using  a  quantum  mechanical  comput¬ 
ing  device  to  efficiently  simulate  quantum  mechanical 
many-body  dynamics  [1-3],  a  task  that  is  exponen¬ 
tially  complex  in  the  number  of  particles  treated  and 
is  completely  intractable  by  any  classical  computing 
means  for  large  systems  of  many  particles.  In  the  two 
decades  following  his  work,  remarkable  progress  has 
been  made  both  theoretically  and  experimentally  in  the 
new  field  of  quantum  computation  [4,5].  Ironically, 
however,  most  of  the  theoretical  progress  in  quan¬ 
tum  computing  has  developed  within  the  purview  of 
the  computer  scientist  with  the  principle  applications 
of  efficient  quantum  information  processing  related  to 
cryptography  and  secure  quantum  communication.*  1  In 
an  effort  return  to  Feynman’s  original  direction,  the 
Air  Force  Research  Laboratory  and  the  Air  Force  Of¬ 
fice  of  Scientific  Research  has  established  a  multidis¬ 
ciplinary  basic  research  theme  called  Quantum  Com¬ 
putation  for  Physical  Modeling  to  explore  quantum  al¬ 
gorithms  to  model  dynamical  physical  systems.  Our 
goal  is  to  establish  a  practical  and  generic  means  by 
which  the  power  of  quantum  mechanics  (that  is,  quan¬ 
tum  parallelism  due  to  the  superposition  and  entangle¬ 
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1  References  to  specific  publications  in  these  subjects  arc  so 

ubiquitous  in  the  quantum  computing  literature  that  we  do  not 
include  any  here.  Comprehensive  collections  of  quantum  computing 
papers  have  been  recently  published  [4,5]. 


ment  of  states)  can  be  used  to  speedup  numerical  sim¬ 
ulations  of  interest  to  computational  physicists. 

Notwithstanding  the  veritable  stampede  towards 
computer  science  related  applications  by  most  re¬ 
searchers  in  the  field  of  quantum  computing,  a  few 
maverick  physicists  have  developed  some  quantum 
algorithms  to  model  quantum  mechanical  systems. 
A  starting  point  for  this  development  was  a  prob¬ 
lem  posed  by  Feynman  himself  to  show  that  the  one 
dimensional  Dirac  equation  could  be  modeled  by  a 
single-speed  particle  traveling  in  a  two-dimensional 
space-time  as  a  sum  over  zigzag  paths  of  straight  line 
elements  [6],  witlyhe  amplitude  of  a  particular  path 
contributing  to  the  kernel  by  the  number  of  “colli¬ 
sions”  or  comers  along  that  zigzag  path.  This  quan¬ 
tum  lattice  gas  representation  of  quantum  mechanics 
is  equivalent  to  the  well  known  path  integral  repre¬ 
sentation.2  A  quantum  lattice  gas  accounts  for  all  con¬ 
tributing  paths  by  simultaneously  evolving  many  par¬ 
ticles  in  a  unitary  fashion.  Therefore,  instead  of  sum¬ 
ming  (or  integrating  over)  paths  as  individual  entities, 
all  contributing  paths  are  effectively  simulated  in  one 
fell  swoop  as  a  combined  field  quantity.  In  the  end,  the 
collisional  interaction  between  particles  in  the  quan¬ 
tum  lattice  representation  can  be  described  by  an  ef¬ 
fective  field  theory  (the  Dirac  equation  in  this  particu¬ 
lar  case)  at  the  large-scale  called  the  continuum  limit. 

Beginning  in  the  mid  1990’s,  a  contemporaneous 
series  of  quantum  lattice-gas  algorithms  to  model 
the  relativistic  Dirac  equation,  equivalent  to  Fcyn- 


2  A  solution  to  Feynman’s  “quantum  lattice  gas”  problem  was 
published  in  1984  by  Jacobson  and  Schulman  [7]. 
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man’s  original  algorithm,  were  published  by  Succi  [8, 
9],  Bialynicki-Birula  [10],  and  Meyer  [11-15].  Fur¬ 
thermore,  a  series  of  papers  on  modeling  the  non- 
relativisitic  Schroedinger  equation  were  published  by 
Boghosian  and  Taylor  [16-18]  and  by  Yepez  and 
Boghosian  [19],  the  latter  article  appearing  in  this  is¬ 
sue.  Our  present  goal  in  the  Quantum  Computation  for 
Physical  Modeling  project  is  accelerate  this  algorith¬ 
mic  developmental  effort  that  has  occurred  over  the 
past  decade. 

In  fact,  we  hope  to  go  further  in  the  application  of 
this  quantum  algorithmic  method.  We  have  developed 
new  efficient  quantum  lattice-gas  algorithms  to  model 
classical  dynamical  systems  [20-23].  Meyer  also  ad¬ 
dresses  this  subject  in  his  article  on  physical  quan¬ 
tum  algorithms  contained  in  this  issue  [24].  In  the 
past,  we  have  considered  quantum  algorithms  suited 
to  globally  phase-coherent  quantum  computers  [25], 
however  our  focus  is  presently  on  those  quantum  al¬ 
gorithms  suited  to  implementation  on  locally  phase- 
coherent  quantum  computers  that  are  technologically 
much  simpler  to  experimentally  implement  [26].  This 
program  to  use  one  quantum  mechanical  system  to 
model  another  quantum  mechanical  system  or  clas¬ 
sical  system  can  perhaps  best  be  described  as  effi¬ 
cient  analog  computing,  which  in  this  context  may  be 
termed  analog  quantum  computing. 

The  principle  technological  obstacle  to  globally 
phase-coherent  quantum  computation  is  the  problem 
of  the  uncontrolled  decoherence  of  the  quantum  com¬ 
puter’s  wave  function.  The  quantum  computing  com¬ 
munity  at  large,  following  the  traditional  computer  sci¬ 
entist’s  mindset  for  correcting  bit-flip  errors  using  re¬ 
dundancy,  has  been  investing  much  theoretical  work 
in  attempts  to  develop  generalized  methods  for  quan¬ 
tum  error  correction  of  bit-flip  and  phase-change  er¬ 
rors  [27-29].  As  an  expedient  alternative  to  this  cum¬ 
bersome  approach,  as  demonstrated  by  the  advent  of 
several  nuclear  magnetic  resonance  quantum  comput¬ 
ing  experiments  [30-33],  it  is  possible  to  avoid  bit-flip 
and  phase-change  errors  altogether:  Keep  the  individ¬ 
ual  quantum  computing  elements  small  enough  so  that 
all  computation  occurs  within  a  single  spin-spin  deco¬ 
herence  time  where  bit-flip  and  phase-change  errors 
are  irrelevant. 

Given  this  possibility  of  avoiding  errors,  it  is  nat¬ 
ural  to  consider  building  a  large-scale  quantum  com¬ 
puting  system  by  connecting  many  small  ones  together 


in  an  array  interconnected  by  nearest-neighbor  clas¬ 
sical  communication  channels.  We  call  this  type  of 
quantum  mechanical  device  a  type-II  quantum  com¬ 
puter  [26].  This  hybrid  architecture,  combining  clas¬ 
sical  massively  parallelism  and  quantum  parallelism, 
is  suited  to  modeling  dynamical  physical  systems, 
such  as  turbulent  Navier-Stokes  fluids  [21-23,25].  In 
collaboration  with  the  Nuclear  Engineering  Depart¬ 
ment  of  MIT,  we  have  developed  a  prototype  type-II 
quantum  computer  based  on  spatial  nuclear  magnetic 
resonance  spectroscopy.  We  use  a  gradient  magnetic 
field  to  distinguish  individual  layers  in  a  liquid  sam¬ 
ple  so  that  each  layer  effectively  becomes  an  indi¬ 
vidual  quantum  computing  node  comprising  an  en¬ 
semble  of  molecules.  The  first  simulation  of  a  quan¬ 
tum  lattice-gas  model  for  the  one-dimensional  diffu¬ 
sion  equation  [22]  has  been  carried  out  on  this  quan¬ 
tum  computer  prototype  using  the  atomic  spin-state  of 
Carbon- 13  and  Hydrogen  nuclei  within  a  linear  array 
of  chloroform  molecules  [34].  This  milestone  repre¬ 
sents  the  first  physical  simulation  to  date  on  a  quantum 
computer  and  is  contained  in  this  issue.  A  subsequent 
paper  presenting  an  improved  version  of  our  type-II 
quantum  computer  prototype,  that  corrects  for  various 
implementation  errors  and  uses  better  quantum  con¬ 
trol,  is  also  in  preparation  [35]. 

The  rather  rapid  proof-of-concept  achieved  by  spa¬ 
tial  nuclear  magnetic  resonance  spectroscopy  of  a 
type-II  quantum  computer  has  led  to  interest  in  the 
subject  by  the  Office  of  the  Secretary  of  the  Air  Force 
and  the  House  Appropriations  Committee  of  the  US 
Congress  leading  in  turn  to  a  strong  commitment  to 
the  field  of  quantum  computation  for  physical  mod¬ 
eling  by  the  Department  of  Defense.  Our  new  quan¬ 
tum  computing  research  theme  is  supported  by  sev¬ 
eral  directorates  of  the  Air  Force  Office  of  Scientific 
Research  with  about  two  dozen  university  research 
projects  across  the  country  to  date.  The  design  and 
construction  of  several  new  type-II  quantum  computer 
prototypes  are  now  underway  using  various  technolo¬ 
gies  including  superconducting  electronics  and  quan¬ 
tum  optics  for  example. 

We  have  established  a  new  annual  workshop  se¬ 
ries  dedicated  to  this  subject  of  quantum  compu¬ 
tation  for  physical  modeling  (see  our  web  site  at 
http://qubit.plh.af.mil  for  more  details).  The  first  work¬ 
shop  in  this  series  was  held  in  the  fall  of  2000  in  North 
Falmouth  in  Cape  Cod,  Massachusetts  and  the  fol- 
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lowing  collections  of  articles  contained  in  this  issue 
were  contributed  from  this  workshop.  Our  goal  for  this 
workshop  series  is  to  annually  publish  a  collection  of 
such  contributed  articles,  to  review  our  progress,  and 
to  provide  an  open  forum  for  new  collaborators  to  join 
us  in  this  activity. 
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Abstract 

Presented  is  quantum  lattice-gas  model  for  simulating  the  time-dependent  evolution  of  a  many-body  quantum  mechanical 
system  of  particles  governed  by  the  non-relativistic  Schrodinger  wave  equation  with  an  external  scalar  potential.  A  variety  of 
computational  demonstrations  are  given  where  the  numerical  predictions  are  compared  with  exact  analytical  solutions.  In  all 
cases,  the  model  results  accurately  agree  with  the  analytical  predictions  and  we  show  that  the  model’s  error  is  second  order 
in  the  temporal  discretization  and  fourth  order  in  the  spatial  discretization.  The  difficult  problem  of  simulating  a  system  of 
fermionic  particles  is  also  treated  and  a  general  computational  formulation  of  this  problem  is  given.  For  pedagogical  purposes, 
the  two-particle  case  is  presented  and  the  numerical  dispersion  of  the  simulated  wave  packets  is  compared  with  the  analytical 
solutions.  ©  2002  Published  by  Elsevier  Science  B.V. 

Keywords:  Schrodinger  wave  equation;  Quantum  computing;  Quantum  lattice  gas;  Quantum  mechanics;  Computational  physics 


1.  Introduction 

Feynman’s  work  regarding  quantum  mechanical 
computers  used  to  simulate  physical  quantum  me¬ 
chanical  behavior  in  a  numerically  efficient  way  [1-3] 
has  subsequently  led  to  several  series  of  research  pa¬ 
pers  concerned  with  a  variety  of  details  involving  the 
particular  quantum  algorithm  that  best  represents  the 
Feynman  path  integral  [4,5].  The  quantum  algorithmic 
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approach  is  based  on  a  two-component  complex  field 
defined  on  a  discrete  spacetime  lattice  where  unitary 
matrices  act  locally  on  the  field  causing  its  temporal 
evolution  in  discrete  time  steps.  Using  such  a  spatially 
discrete  field  makes  it  possible  to  computationally  rep¬ 
resent,  in  the  long  wavelength  limit  of  modes  in  the 
discrete  system,  the  dynamical  time-dependent  evolu¬ 
tion  of  a  continuous  wave  function  in  a  manner  that  is 
numerically  efficient. 

In  1 994  Bialynicki-Birula  presented  a  general  quan¬ 
tum  algorithmic  approach  of  this  kind  for  modeling 
the  Weyl,  Dirac,  and  Maxwell  equations  on  a  body- 
centered  cubic  lattice  in  three-dimensions  [6].  In  a 
series  of  papers  on  simulating  the  one-dimensional 
Dirac  equation  [7-9],  Meyer  presented  a  quantum  al¬ 
gorithm  similar  to  that  of  Bialynicki-Birula  with  a  va- 
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riety  of  numerical  simulations  including  the  effects  of 
boundary  conditions,  inhomogeneities,  and  an  exter¬ 
nal  scalar  potential.  Meyer  set  the  quantum  algorithm 
for  the  discretized  path  integral  in  the  context  of  what 
is  called  the  quantum  lattice  gas  method  and  his  algo¬ 
rithm  is  equivalent  to  the  one-dimensional  version  of 
the  Bialynicki-Birula  quantum  algorithm  for  the  Dirac 
equation. 

Contemporaneously  with  Meyer,  two  other  series 
of  papers  on  quantum  lattice-gas  models  of  the  one¬ 
dimensional  non-relativistic  Schrodinger  wave  equa¬ 
tion  were  presented  by  Succi  and  Benzi  [10,11]  and 
by  Boghosian  and  Taylor  [12-14].  The  approach  taken 
by  Succi  and  Benzi  is  somewhat  more  computation¬ 
ally  oriented  in  that  they  begin  with  a  “kinetic”  lat¬ 
tice  Boltzmann  equation  of  motion1  (effectively  the 
one-dimensional  Dirac  equation  in  the  Majorana  rep¬ 
resentation)  and  show  that  the  Schrodinger  wave  equa¬ 
tion  emerges  as  the  governing  equation  of  motion  for 
the  slow  mode  in  the  long  wavelength  “hydrodynam¬ 
ic”  limit.  That  is,  Succi  and  Benzi  observed  that  the 
“macroscopic  scale”  Schrodinger  wave  equation  arises 
from  the  “mesoscopic  scale”  Dirac  equation  in  a  man¬ 
ner  quite  analogous  to  how  the  macroscopic  Navier- 
Stokes  hydrodynamic  fluid  equation  arises  from  the 
mesoscopic  kinetic  Boltzmann  equation  through  the 
Chapman-Enskog  expansion. 

Boghosian  and  Taylor’s  approach  follows  more 
along  the  lines  of  Meyer’s  approach  in  that  their  model 
is  developed  as  a  generalization  of  the  classical  lattice 
gas  method.  A  kinetic  transport  equation,  now  formu¬ 
lated  directly  at  the  “microscopic  scale”,  again  leads 
to  the  Schrodinger  wave  equation  in  the  continuum 
limit.  The  Boghosian  and  Taylor  quantum  lattice  gas 
model  focuses  on  solving  the  many-body  Schrodinger 
wave  equation  with  an  arbitrary  scalar  potential  in  an 
arbitrary  number  of  spatial  dimensions.  They  analyti¬ 
cally  argue  for  an  exponential  numerical  speedup  aris¬ 
ing  from  simulation  in  the  many-body  sector  of  the  full 
Hilbert  space  carried  out  simultaneously  using  quan¬ 
tum  superposition  of  states.  The  Boghosian  and  Taylor 
version  of  the  quantum  algorithm  is  also  cast  explic¬ 
itly  for  direct  implementation  of  an  array  of  quantum 
bits  [13].  Polley  has  recently  presented  an  argument 


1  The  classical  lattice  Boltzmann  equation  was  popularized  with 
its  application  to  computational  fluid  dynamics  [15-17]. 


for  adding  both  an  external  scalar  and  vector  poten¬ 
tial  into  a  quantum  lattice-gas  model  by  analytically 
demonstrating  the  discrete  model’s  invariance  with  re¬ 
spect  to  a  general  local  gauge  transformation  [18]. 

A  characteristic  feature  of  all  these  quantum  algo¬ 
rithms,  used  to  model  the  dynamical  behavior  of  ei¬ 
ther  relativistic  or  non-relativistic  quantum  particles, 
is  that  the  governing  wave  function  is  well  approxi¬ 
mated  as  one  approaches  the  continuum  limit  where 
the  grid  resolution  of  the  spatial  lattice  become  in¬ 
finite  (the  lattice  cell  size  approaches  zero).  There¬ 
fore,  from  the  point-of-view  of  the  modeler,  there  ex¬ 
ists  a  “microscopic  scale”  where  the  unitary  dynamical 
rules  are  locally  applied  in  a  discretized  fashion  and 
time  advances  forward  in  incremental  units  in  a  way 
that  is  quite  artificial.  Yet  at  the  “macroscopic  scale”, 
which  corresponds  to  the  long  wavelength  limit  of  the 
dynamical  modes  in  the  discrete  system,  there  is  an 
emergent  effective  field  theory  for  a  complex  ampli¬ 
tude  field,  continuous  and  differentiable  in  both  space 
and  time,  that  exactly  obeys  the  physical  quantum  me¬ 
chanical  equations  of  motion.  At  the  small  scale  (char¬ 
acterized  by  the  lattice  cell  size)  one  imagines  a  sys¬ 
tem  of  fermionic  particles  undergoing  local  collisions 
and  translation  to  nearby  nodes  of  the  lattice.  Each  of 
these  fermionic  particles  occupies  a  local  positional 
state  at  a  specific  lattice  node  with  a  certain  probabil¬ 
ity  amplitude.  All  the  possible  locations  of  the  actual 
physical  quantum  particle  are  effectively  modeled  by 
the  interfering  set  of  probability  amplitudes  associated 
with  this  system  of  fermionic  particles.  That  is,  all  the 
possible  pathways  are  modeled  simultaneously  using  a 
kind  of  kinetic  system  of  locally  interacting  fermions 
on  the  small  scale. 

Seen  as  a  kinetic  system  then,  we  may  expect  that 
there  exists  a  local  equilibrium  configuration  of  parti¬ 
cles.  We  require  that  this  configuration  be  an  eigen- 
ket,  with  unity  eigenvalue,  of  the  local  unitary  col¬ 
lision  operator  that  is  uniformly  and  spatially  homo¬ 
geneously  applied  to  the  lattice-based  two-component 
field.  Then  the  dynamical  lattice-based  system  on  all 
lattice  nodes  undergoes  local  relaxation  towards  this 
equilibrium  configuration.  However,  unlike  a  classi¬ 
cal  kinetic  system,  the  global  configuration  of  particles 
does  not  relax  to  a  single  steady-state  equilibrium.  In¬ 
stead,  there  arc  many  steady-state  global  equilibrium 
configurations  which  effectively  are  the  energy  eigen¬ 
states  of  the  “macroscopic  scale”  quantum  mechani- 
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cal  equation  of  motion  modeled  by  the  quantum  lattice 
gas.  If  the  quantum  lattice  gas  system  is  initialized  in 
any  one  of  the  energy  eigenstate  global  configurations, 
the  macroscopic  scale  configuration  of  the  system  will 
remain  fixed  in  time  albeit  the  microscopic  configu¬ 
ration  of  particles  would  be  continually  changing  at 
every  unit  time  step.2  In  the  end,  the  Feynman  path  in¬ 
tegral  is  efficiently  and  accurately  recast  as  a  kinetic 
dynamical  process  computed  in  parallel  efficiently  on 
a  spacetime  lattice. 

In  this  paper  we  do  not  argue  that  the  quantum 
lattice  gas  dynamical  rules  represent  a  local  hidden 
variable  theory  of  quantum  mechanics.  Although  fun¬ 
damental  arguments  can  be  made  to  limit  the  possi¬ 
ble  form  of  the  local  unitary  operator  in  the  quan¬ 
tum  lattice  gas  [6,7],  these  arguments  lead  to  an  al¬ 
gorithm  suited  for  implementation  on  a  quantum  com¬ 
puter.  In  our  investigation  of  a  suitable  local  unitary 
collision  operator  we  have  found  that  one  quantum 
gate  in  particular,  the  square-root-of-swap  given  be¬ 
low  in  Section  2.3,  is  especially  useful  for  model¬ 
ing  the  Schrodinger  wave  equation  in  that  the  local 
equilibrium  configuration  discussed  above  is  an  eigen- 
ket  of  this  gate  and  with  unity  eigenvalue.  There¬ 
fore,  we  have  selected  the  square-root-of-swap  gate 
as  our  basic  model  quantum  gate.  As  demonstrated  in 
Section  2.6,  we  find  that  this  quantum  gate  leads  to 
an  overall  modeling  error  that  is  second  order  in  the 
temporal  discretization  and  fourth  order  in  the  spatial 
discretization.  Another  point  regarding  this  particular 
gate  is  that  when  measurements  are  periodically  made 
of  the  state  of  qubits  in  the  system,  which  destroys 
quantum  superpositions  and  entanglements  in  the  sys¬ 
tem,  the  macroscopic  scale  behavior  of  the  quantum 
lattice  gas  system  is  governed  by  the  classical  diffu¬ 
sion  equation  [19]. 

This  paper  is  organized  as  follows.  First,  in  Sec¬ 
tion  2  we  describe  the  basic  quantum  algorithm,  in 
particular,  how  we  encode  the  wave  function,  how  we 
use  two  basic  quantum  gates  applied  to  the  qubits  at 
the  nodes  of  the  lattice.  We  first  describe  the  algo¬ 
rithm  using  matrices  and  then  we  describe  an  equiv¬ 


2  Given  a  finite  size  lattice  used  for  modeling  purposes,  each  lo¬ 
cal  configuration  oscillates  in  time  even  when  the  global  configura¬ 
tion  is  a  time-independent  energy  eigenvalue.  However,  the  ampli¬ 
tude  of  the  oscillation  does  approach  zero  as  the  lattice  size  becomes 
infinite. 


alent  finite  difference  formulation.  Using  the  micro¬ 
scopic  finite  difference  equations,  we  then  derive  the 
effective  field  theory  at  the  macroscopic  scale  where 
both  the  lattice  cell  size  and  the  update  time  step  ap¬ 
proach  zero.  Diffusive  ordering  holds,  as  is  typical  for 
lattice  gas  systems,  where  small-scale  temporal  fluc¬ 
tuations  in  the  wave  function  go  as  the  square  of  the 
magnitude  of  the  small-scale  spatial  fluctuations.  To 
confirm  that  our  derivation  is  correct  and  to  test  the 
validity  of  the  quantum  lattice  gas  model,  we  test  the 
time-dependent  dispersion  of  a  free  Gaussian  packet. 
We  also  test  the  system  when  it  is  initialized  in  an  en¬ 
ergy  state,  which  is  a  fixed  macroscopic  scale  steady- 
state  configuration.  We  find  that  as  we  halve  the  lattice 
cell  size  (double  the  spatial  resolution)  the  cumulative 
error  in  the  model  drops  by  a  factor  of  25  45  =  43.7. 

Second,  in  Section  3  we  show  how  an  external 
scalar  potential  may  be  modeled  in  the  quantum  lattice 
gas  system  by  inducing  a  local  phase  rotation  in  the 
qubits  at  each  node  of  the  lattice.  The  qubits  at  a  lattice 
node  are  phase  rotated  by  an  amount  corresponding  to 
the  strength  of  the  spatially  dependent  external  poten¬ 
tial  at  that  node.  We  show  how  this  local  phase  rota¬ 
tion,  a  kind  of  gauge  transformation,  produces  an  ad¬ 
ditional  potential  energy  term  in  the  Schrodinger  wave 
equation.  We  then  test  the  quantum  algorithm  against 
two  well-known  cases  of  harmonic  oscillation  in  a  par¬ 
abolic  potential  well  and  quantum  scattering  off  of  and 
tunneling  through  a  constant  potential  energy  barrier. 
In  both  cases,  the  model  behaves  as  expected. 

Third,  in  Section  4  we  test  how  well  the  quantum 
lattice  gas  can  simulate  the  simultaneous  dispersion  of 
two  fermionic  particles.  In  the  case  of  multiple  parti¬ 
cles,  the  operational  quantum  gate  sequence  to  handle 
the  many-body  case  is  identical  to  the  single  particle 
case  presented  in  Section  2.  Therefore,  the  implemen¬ 
tation  of  cither  situation  on  a  quantum  computer  would 
be  identical  as  well,  except  for  state-preparation  and 
measurements.  However,  since  at  present  no  quantum 
computer  exists  that  can  test  the  algorithm  presented 
here,  we  are  forced  to  consider  the  implementation 
on  a  classical  computer  and  here  the  implementation 
for  a  many-body  problem  is  much  more  complex  than 
for  the  single-body  problem.  Nevertheless,  we  present 
a  general  formulation  of  the  quantum  gates  in  a  sec¬ 
ond  quantized  representation  where  the  basic  compu¬ 
tational  operations  are  creation  and  annihilation  of  lo¬ 
cal  particle  occupancies.  The  advantage  of  such  an 
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implementation  is  that  it  is  straightforward  to  imple¬ 
ment  the  fundamental  creation  and  annihilation  opera¬ 
tion  is  a  way  that  respects  the  anti-commutation  rela¬ 
tions  for  fermionic  particles.  We  demonstrate  that  the 
macroscopic  scale  behavior  of  the  quantum  lattice  gas 
agrees  with  the  exact  time-dependent  solution  of  the 
two-body  wave  equation. 

Finally,  we  conclude  this  paper  with  a  short  de¬ 
scription  of  results  and  share  some  lessons  we  learned 
after  using  the  quantum  lattice  gas  model  extensively. 
We  also  point  out  some  future  directions  that  may  be 
taken  to  expand  the  usefulness  of  the  model. 

2.  Quantum  algorithm  for  a  single  free  particle 

We  describe  the  quantum  lattice-gas  algorithm  for 
modeling  the  Schrodinger  wave  equation  by  consid¬ 
ering  the  simplest  case  of  a  single  free  particle  in  a 
one-dimensional  space.  In  this  simple  case,  the  wave 
function  /)  obeys  the  following  partial  differen¬ 
tial  equation  in  the  position  representation 

.£»*(*. o  h2  a2*(*.o  /n 

'  dt  2m  dx2  ’ 

where  h  is  Planck’s  constant  and  m  is  the  mass  of 

the  quantum  particle.  Here  \//(x ,  f)  is  a  continuous 

probability  amplitude  field  (e.g.,  a  continuous  complex 

field). 

2. 1.  Encoding  the  wave  function 

To  “program”  a  quantum  computer  to  simulate  (1), 
it  is  necessary  to  first  formulate  an  encoding  scheme 
where  a  collection  of  qubits  is  used  to  store  the  value 
of  the  wave  function.  Since  the  number  of  qubits  in 
any  quantum  computer  is  necessarily  a  finite  number, 
the  wave  function  will  have  to  be  approximated  in 
the  usual  way  by  discretizing  a  physically  continuous 
amplitude  field  into  an  artificially  discrete  and  finite 
set  of  complex  numbers.  To  do  this,  let  us  begin  with 
a  one-dimensional  spatial  lattice  with  L  number  of 
nodes.  With  each  node  of  the  lattice  we  associate  a 
position  basis  ket  denoted  by  |jc/>,  where  0  ^  /  ^ 
L  —  1 .  The  discretized  system  ket  in  the  position  basis 
is 

L- 1 

|  yo  =  £>!*/)■  (2) 

1=0 


where  c/  =  (jc/ \\j/)  is  a  complex  number.  In  other 
words,  the  basic  approach  to  model  the  single  particle 
wave  function  governed  by  (1)  is  to  express  \x[/)  as 
a  sum  of  all  the  possible  ways  the  particle  can  be 
situated  on  the  lattice  with  a  probability  amplitude  c/ 
associated  with  each  possible  location  \x /). 

In  our  model,  we  assign  two  qubits  to  each  node 
of  the  lattice,  for  a  total  of  2 L  qubits  in  the  whole 
quantum  computer.  The  qubits  that  reside  at  the  / th 
node  of  the  lattice  are  denoted  by  \ql0)  and  \q\)  and 
they  are  used  to  encode  the  coefficient  c/  of  (2)  of 
the  position  ket  for  that  node.  Each  qubit  is  a  two- 
level  quantum  system  \qla)  =  otla |0)  -I-  Pla\l),  where 
\otla\ 2  4-  \fla\2  =  1  for  a  =  0  or  1  and  0  ^  ^  L  -  1. 
We  consider  each  qubit  to  be  a  container  that  may 
or  may  not  be  occupied  by  the  quantum  particle.  The 
quantum  particle  is  said  to  occupy  the  #  th  local  state  at 
position  xi  when  /3la  =  1.  Similarly,  the  a  th  local  state 
at  position  x /  is  said  to  be  empty  when  =  0. 

To  see  how  the  qubit  encoding  works,  we  write 
\yjr)  in  the  number  representation.  In  the  number 
representation,  each  basis  state  is  expressible  as  the  ket 
‘  ‘  'non\  )>  where  nla  =  0  or  1  for  all  / 
and  a.  The  Boolean  variables  nla  are  called  the  number 
variables  and  they  correspond  to  a  binary  indexing  of 
the  basis  states  in  the  number  representation.  Since  we 
are  concerned  with  modeling  the  one-particle  wave 
equation,  we  need  consider  only  a  subset  of  all  the 
basis  states  where  only  one  of  the  number  variables 
is  1  and  all  the  rest  are  0.  This  subset  of  all  the  basis 
states  is  called  the  one-particle  sector.  There  are  2 L 
such  combinations  and  we  shall  label  these  with  the 
binary  encoding  formula  |22/_N/),  for  a  =  0,  1  and 
0  ^  ^  L  —  1 .  Therefore,  the  system  ket  in  the  number 
representation  can  be  written  as 

L- 1  1 

W  =  £!>/+*  |22/+a>,  0) 

/=0  a=  0 

where  each  %21+a  is  a  probability  amplitude  (e.g., 
complex  number). 

Now  for  each  position  ket  |jc/)  there  are  two  cor¬ 
responding  basis  states  in  the  number  representation 
|22/)  and  |22/+1).  There  are  two  interfering  possibil¬ 
ities  for  a  particle  to  occupy  the  /th  position  on  the 
lattice.  Therefore,  the  occupancy  probability  of  the  /th 
node  is  computed  by  first  summing  the  probability  am¬ 
plitudes  of  these  corresponding  basis  states  and  then 
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computing  the  square  of  the  absolute  value  thereof.  In 
other  words,  the  coefficient  c/  in  (2)  is  set  equal  to  the 
sum  of  the  on-site  coefficients  in  (3) 

Q  =£2/+ &/+1-  (4) 


damard  gate  at  the  very  end  of  the  simulation  prior 
to  making  a  measurement  of  the  wave  function  so 
that  a  single  qubit  at  each  node  will  encode  the 
probability  amplitude  c/  =  (jc/  |^r)  in  (2). 


The  definition  (4)  is  an  essential  part  of  the  quantum 
lattice-gas  model  presented  in  this  paper.  In  the  section 
below,  where  we  analytically  predict  an  effective 
field  theory  for  our  artificially  discretized  model,  we 
explain  why  we  need  to  make  this  assignment.  We 
will  find  that  (4)  is  needed  for  the  predicted  effective 
field  theory  to  accurately  approximate  the  Schrodinger 
wave  equation  in  the  long-wavelength  limit,  which  is 
also  defined  below. 

2.2.  Formulating  a  suitable  gate  sequence 

We  shall  require  that  the  algorithmic  scheme  be 
at  least  second  order  convergent  in  space,  so  that 
as  we  double  the  grid  resolution  (e.g.,  double  the 
number  of  qubits)  we  in  turn  reduce  the  numerical 
error  due  to  the  field  discretization  by  a  factor  of  one- 
quarter.  With  this  type  of  convergence  characteristic, 
we  are  assured  that  we  can  simulate  a  wave  function 
governed  by  the  Schrodinger  wave  equation  (1)  to  any 
arbitrary  degree  of  accuracy.  After  we  formulate  our 
algorithmic  scheme,  we  will  then  a  posteriori  verify  by 
direct  numerical  simulation  that  it  is  indeed  at  least  a 
second-order  convergent  numerical  scheme.  In  fact  in 
Section  2.6  we  will  find  that  our  numerical  scheme  is 
fourth-order  convergent  with  an  error  that  goes  as  8x1 2 3  4. 

To  simulate  the  quantum  behavior  of  the  wave 
function,  we  seek  to  develop  a  sequence  of  2-qubit 
gate  operations  that  will  act  on  a  large  collection  of 
qubits  in  the  simplest  way.  We  impose  the  following 
four  simplifying  constraints: 

(1)  All  quantum  gate  operations  are  homogeneous 
and  independent  of  space  and  time. 

(2)  Only  a  single  quantum  gate  is  used  to  evolve  the 
wave  function  and  this  gate  is  applied  to  each 
lattice  node  independently  (locality). 

(3)  To  provide  communication  channels  between  lat¬ 
tice  nodes,  only  the  simplest  gate  is  used  (e.g.,  a 
swap  gate). 

(4)  Because  the  final  value  of  the  computed  wave 
function  depends  on  summing  interfering  possi¬ 
bilities  according  to  (4),  we  shall  use  the  Had- 


With  two  qubits  per  node,  there  are  four  on-site  ba¬ 
sis  kets,  |O)0|O)  =  (1,0,  0,0),  |O)0|1)  =  (0,  1,0,0), 
|1)  <g>  |0)  =  (0,0,  1,0),  and  1 1)  0  1 1)  =  (0,  0,  0,  1).  In 
the  context  of  a  quantum  lattice-gas  model,  the  unitary 
matrix  U  is  called  the  local  collision  operator  and  the 
on-site  ket  \  v)  =  |0)  0  1 1)  +  1 1)  0  |0)  =  (0,  1 ,  1 , 0)  is 
called  the  number  density  ket.  To  have  a  well  behaved 
local  equilibrium  associated  with  the  collision  process, 
the  local  collision  operator  must  have  the  number  den¬ 
sity  ket  as  an  eigenvector  with  unity  eigenvalue. 

2.3.  Matrix  representation 


The  quantum  gate  that  we  use  to  evolve  the  wave 
function,  which  is  applied  independently  on  a  site-by¬ 
site  basis,  is  the  square-root-of-swap  gate 


/I  0  0  0  \ 

0  3  “3  3  +  3  0 

0  3  +  2  3  —  3  0 

VO  0  0  -1/ 


(5) 


The  reason  for  calling  this  the  square-root-of-swap 
gate  is  that  U  2  is  the  swap  gate  itself 


U  U  = 


/I  0  0  0\ 
0  0  10  1 
0  10  0 
\0  0  0  1/ 


(6) 


The  two  nontrivial  eigenvalues  of  U  are  k\  =  1 
and  A 2  =  — i,  with  eigenvectors  |vi)  =  (0,  1,  1,0) 
and  |  V2>  =  (0,  —  1,  1, 0),  respectively.  Also,  since  (5) 
causes  mixing  only  between  the  single-particle  basis 
kets  |0)  0  |1)  and  |1)  0  |0),  it  conserves  particle 
number.  So  (5)  is  an  appropriate  choice  for  the  local 
collision  operator. 

The  full  collision  operator,  denoted  C,  which  acts 
on  the  system  ket  \yfr)  is  formed  by  a  L-fold  tensor 
product  over  the  local  collision  operators  U  applied 
homogeneously  and  independently  on  each  node  of 
the  lattice 


L- 1 

c  =  0Z7. 


1=0 


(7) 
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Let  us  denote  the  swap  operator  by  x/XtV,  where 
pi  and  v  index  any  two  qubits  in  the  system.  The 
streaming  operator,  denoted  Si,  causes  a  global  shift 
to  the  right  of  the  first  qubit  on  all  the  lattice  nodes. 
Therefore,  Si  can  be  represented  by  a  sequence  of 
swap  operators  acting  on  nearest  neighbors 

(L- 1)/2 

S\  =  ]~~[  X2/.2/+2-  (8) 

1=0 

In  matrix  form,  x  is  a  22  x  22  permutation  matrix 
/I  0  0  0  \ 

X  =  °  °  1  °  .  (9) 

*0100  w 

Vo  0  0  -1/ 

The  algorithm  we  use  to  model  the  Schrodinger  wave 
equation  involves  multiple  applications  of  the  collision 
operator  interleaved  with  streaming  operations  as 
follows: 

\t(t  +  r/2))  =  El/2\t(t)),  (10) 

where  the  square  root  of  the  evolution  operator  is 
£,i/2  =  5,tC5iC.  (11) 

Here  S,T  denotes  the  transpose  of  Si  and  is  the  inverse 
of  Si .  Application  of  S,T  causes  a  global  shift  to  the 
left  of  the  first  qubit  on  all  the  lattice  nodes.  One  full 
time  step  of  the  evolution  is 

\r/,(t  +  T))  =  El\x/r(0).  (12) 

We  use  four  applications  of  the  collision  operator  in 
E\  because  C4  is  the  identity  operation.  Note  that  Si 
and  C  do  not  commute,  otherwise  (12)  would  be  a 
trivial  evolution  equation. 

Note  that  in  (12),  our  choice  of  streaming  the  first 
qubit  was  arbitrary.  A  more  balanced  algorithmic  ap¬ 
proach  would  treat  both  qubits  identically.  Therefore, 
we  could  alternatively  define  one  full  time  step  as 

\lr(t  +  T))  =  E2E\\*(t)),  (13) 

where 

e}'2  =  S?CS2C,  (14) 

and  where  the  streaming  operator  S2  causes  a  global 
shift  to  the  right  of  the  second  qubit  on  all  the  lat¬ 
tice  nodes.  The  advantage  of  using  the  balanced  al¬ 
gorithm  (13)  is  that  its  error  is  fourth-order  in  space 
whereas  for  the  unbalanced  algorithm  (12)  it  is  only 
third-order. 


2.4.  Finite  difference  formulation 

It  is  possible  to  specify  the  quantum  algorithm  to 
model  the  Schrodinger  equation  without  the  use  of 
matrices.  Instead  we  can  write  down  a  set  of  finite 
difference  equations,  which  are  equivalent  to  (10), 
but  perhaps  simpler  to  comprehend  at  first  glance. 
To  do  this,  let  us  introduce  a  new  notation  for  the 
2 L  probabilities  amplitudes  £2 i+a  in  (3).  We  will 
denote  the  two  complex  numbers  per  lattice  node  by 
cpo(xi ,  tn)  and  <pi  (*/,  f„).  That  is,  we  have  L-pairs 
of  complex  numbers.  Then,  the  quantum  algorithmic 
operations  (14)  can  be  expressed  as  follows: 

if  mod(/i,  4)  =  0, 

(po(xi,tn)  =  A*<po(xi,  tn-\)  +  A<p\(xi ,  /„_i), 

(15) 

<P\(xt,tn)  =  A<po(x/,  tn-\)  +  A*<p\(xi,  t„- 1), 
if  mod(n,  4)  =  1, 

<P0(xi,tn)  =<poC*/-l,*n-l), 

(16) 

<P\(Xh*n)  =<P\(xi>tn-\), 
if  mod(>i,  4)  =  2, 

(PoixiJn)  =  A*<po(xi,tn-\)  +  A<p\(xi,  tn-\)9 

„  (17) 

<p\(xt,tn)  =  Acpo(xt ,  tn- 1)  -f  A  < p\(xi ,  f„_i), 
and  if  mod(«,  4)  =  3, 

<Po(xi,tn)  =  <PoU/+i,6i-i), 

(18) 

<P\(Xl,tn)=<P\(Xl,tn-\)9 

where  A  =  |  -f  j.  The  finite-difference  equation 
pair  (15)  is  equivalent  to  the  local  collision  operation 
C,  as  is  the  pair  (18).  The  equation  pairs  (16)  and  (18) 
are  equivalent  to  the  streaming  operations  S  and  ST, 
respectively.3 

This  finite-difference  representation  of  the  algo¬ 
rithm  is  nearly  identical  to  that  presented  by  Boghosian 
and  Taylor  in  1997  [14]  where  the  two  on-site  qubits 

3  Noting  that  A  +  A*  =  \  ,  this  set  of  finite  difference  equations 
can  be  expressed  in  a  more  compact  way 

<P0(xi+€Jn+\)=<P0(xhtn)  + 

<P\(xiJn+\)=(p\(xiJn)  +  n\, 

where  €  =  (—1)”  and  Qq  =  A({p\  —  vq)  and  Q\  =  — which  has 
the  standard  form  of  a  lattice-gas  transport  equation. 
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are  simultaneously  streamed  to  the  left  and  right  after 
collision  operation 

<Po(*/+l.f#i)  =  A* <po(xi,  tn-\)  -  A(p\(xi ,  tn- 1), 

(19) 

<P  1  (*/-l»  tn)  =  -A<po(xi,  tn-\)  +  (*/, 

They  noted  that  after  four  time  steps,  the  total  am¬ 
plitude  \lr(xi,tn)  =  (f>o(xi,  tn)  4-  <t>\  (xiJn)  satisfies  a 
finite-difference  equation  which  approximates  the 
Schrodinger  equation  in  the  continuum  limit.  The 
two  essential  differences  between  the  improved  algo¬ 
rithm  (15)  through  (18)  presented  in  this  paper  and 
the  quantum  algorithm  (19)  appearing  in  [20]  is  that 
we  have  alleviated  the  problem  of  the  occurrence  of 
two  non-interpenetrating  lattice-gas  systems  indepen¬ 
dently  evolving  on  different  checker-board  sub-lattices 
and  we  have  doubled  the  numerical  accuracy.  This  is 
a  problem  that  occurs  when  both  on-site  qubits  are  si¬ 
multaneously  streamed  because  streaming  only  a  sin¬ 
gle  qubit  at  a  time,  as  was  done  for  the  quantum  lattice- 
gas  model  of  the  diffusion  equation  [19],  causes  inter¬ 
actions  between  all  the  qubits  at  each  time  step. 

2. 5.  The  governing  partial  differential  equation 

It  is  straightforward  using  a  symbolic  mathematics 
program,  and  tedious  by  hand,  to  use  the  update  rules 
(15)  through  (18)  to  algebraically  determine  the  value 
of  <po  and  tp\  at  a  later  time.  With  the  initial  wave 
function  set  at  to,  one  complete  cycle  of  the  algorithm 
is  completed  at  t%  (that  is,  t%  —  to  =  8t).  With  the  wave 
function  defined  as  ^r(xt,tn)  =  <po(xt,  tn)  +  <p\(xi,tn), 
the  result  after  one  cycle  is4 

'l'(xiJs) 

1  ~b  i 

= - 2~i^(xi,t o)  4-  lM*/4l.  to)  +  o) 

-  [^(*/+2.  to)  +  to)]-  (20) 

Note  that  (20)  is  the  simplified  form  of  the  finite- 
difference  equation  at  the  macroscopic  scale  when  the 
system  is  very  close  to  local  equilibrium  throughout 


the  course  of  the  evolution  as  <po(x,t)  =  <p\(x,t)  ~ 
t)  for  all  x.  The  full  finite-difference  equation 
is  too  long  to  present  here,  but  is  given  in  Appendix  A. 
This  result  is  a  finite-difference  equation  for  the 
following  partial  differential  equation  governing  the 
continuous  amplitude  field  \//(x ,  r) 


d\//(x ,  t) 
dt 


1  8x2  d2\/f(x,  t) 

)  = - 

2  it  dx2 


+  0(8  x4), 

(21) 


which  is  an  approximation  of  (1)  where  the  diffusion 
constant  is  h/m  =  8x2/8t  and  where  8x  is  the  lattice 
cell  size. 

If  one  adds  a  phase  angle  £  to  the  off-diagonal 
components  of  collision  operator  (5)  to  obtain  a 
slightly  more  general  collision  operator 


/I  0  0  0  \ 

0  i-j  (i  +  i)e'<  0 

<•  (5  +  i)'-,(  i~i  0 

Vo  o  o  -1/ 


(22) 


then  the  resulting  governing  partial  differential  equa¬ 
tion  will  have  its  transport  coefficient  dependent  on 
this  phase  angle  as  follows: 


djr(x,t) 

dt 


+  0(8t2) 

i  8x2  d2\j/(x,  t) 
2secf  8t  dx2 


+  0(8x 4). 


(23) 


This  allows  us  to  simulate  a  quantum  system  where 
a  particle’s  mass  can  be  arbitrarily  large  m  =  secf. 
Note  that  in  this  case  the  error  is  cubic  and  is 
proportional  to  sin£.  So  for  very  large  masses,  the 
accuracy  of  the  model  is  reduced  to  third-order  in 
space.  Note  that  (23)  is  is  valid  effectively  field  theory 
at  the  macroscopic  scale  when  the  system  is  very 
close  to  local  equilibrium  where  $po(jc,  t)  =  cp\  (jc,  t) 
^t/c(jc,  t)  for  all  jc. 


2. 6.  Numerical  confirmations 


4  Note  that  the  result  (20)  is  accurate  up  to  fourth  order  in  8x 
only  in  the  situation  where  the  initial  system  is  in  local  equilibrium 
defined  by  (fo(xiJn)  =  ip\{*lJn)-  In  the  more  general  situation 
when  the  system  is  not  in  local  equilibrium  where  <po(xi,tn)  ^ 
<P\  (•*/ .  in  h  the  result  (20)  is  accurate  only  up  to  third  order  in  fix. 


To  numerically  test  that  the  quantum  algorithm  ( 1 2) 
is  indeed  equivalent  to  the  finite-difference  equa¬ 
tion  (20)  and  to  see  just  how  good  of  an  approxima¬ 
tion  of  the  single-particle  Schrodinger  equation  it  is, 
we  have  performed  two  simulations. 
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In  the  first  simulation,  we  test  the  numerical  time 
evolution  of  a  Gaussian  packet 


W*,0) 


1  -*2/(2  a2) 

cr 1  /27T  I/4  C 


(24) 


where  i  ^  x  ^  L  for  a  lattice  of  size  L  =  64£  and 
where  the  packet  width  is  a  =  L/10  as  shown  in 
Fig.  1. 

The  exact  analytical  solution  of  (21)  is  obtained  by 
computing  the  Fourier  components  of  the  energy  basis 
functions 


ao 


an 


L/2 

=  T  I 

-L/2 

L/2 

=  Z  /  ^^’O)cos(^2'77r^ 

-L/2 


dv. 


(25) 


(26) 


With  h  =  1  and  m  =  1 ,  the  energy  eigenvalues  are 
2 n2n2Sx2 

En  — 


m,  •  ,27) 

and  the  time-dependent  solution  to  (21)  plotted  in 
Fig.  1  is 


10  20  30  40  50  60 

Position 


Fig.  1 .  Time  evolution  of  a  Gaussian  packet  for  a  single  quantum 
particle  overplotted  in  succession  where  the  jc -axis  is  the  position 
on  a  64-node  lattice  in  units  of  the  lattice  spacing  i  and  the  y- axis 
is  the  probability  density  \\f/(x,t)\^.  The  solid  curves  are  the  exact 
analytical  solution  and  the  circles  are  the  data  from  the  quantum 
lattice-gas  simulation  (the  initial  wave  function  was  normalized, 
therefore  the  area  under  each  curve  is  one).  The  lattice  size  is 
L  =  64 1.  The  initial  Gaussian  packet  of  with  a  —  L/10  at  t  =  0 
is  centered  at  x  =  32 i  and  the  dispersion  is  evident  by  observing 
the  wave  function  at  the  subsequent  times  t  =  50r,  lOOr,  150r, 
and  200r.  Periodic  boundary  conditions  were  used  and  nm ax  =  20 
energy  eigenmodes  were  used  to  generate  the  exact  solutions.  A  time 
scale  factor  ts  =  1 .04  was  used  to  improve  the  agreement  between 
the  numerical  and  analytical  solutions. 


"ma*  /  x  \ 
V'exact (•*>  0  =  ao  +  an  COs(2n7T  - 


(28) 


Note  that  in  (28),  time  is  scaled  by  a  factor  ts  to 
account  for  kinetic  corrections  to  the  time  step.  As 
the  number  of  lattice  nodes  becomes  large,  this  scaling 
factor  approaches  one. 

The  second  test  of  the  quantum  lattice-gas  algo¬ 
rithm  as  a  model  of  the  Schrodingcr  wave  equation  is 
the  measurement  of  its  numerical  convergence.  Mul¬ 
tiple  simulations  (10  in  total)  were  carried  out  for 
lattice  sizes  ranging  from  L  =  8f ,  16f,  32 f , ...  up  to 
L  =  8192f.  In  each  case  the  initial  state  of  the  simu¬ 
lation  was  the  ground  state  (a  sinusoidal  energy  eigen¬ 
state) 


\j/(x,0)  =  V^exactU)  = 


cos(2nx/L) 

7Z72 


(29) 


Each  simulation  was  run  for  one  time  step  T  =  z  and 
the  numerical  error,  denoted  6,  from  the  exact  solution 
was  then  measured  using  the  following  formula 


e(L)  =  j-  f:{mX,T)\2-\r™'W\2}2.  (30) 

L\J,=. 

We  define  the  grid  resolution  as  the  inverse  of  the  total 
number  of  lattice  points.  That  is,  for  a  box  of  size  1 , 
the  resolving  cell  size  is  defined  as  8x  =  j.  A  plot 
of  the  error  versus  the  resolution  is  given  in  Fig.  2. 
As  the  resolution  is  increased,  the  error  drops  off  as 
€(L)~Z,-5  45. 


Fig.  2.  Log-log  plot  of  the  numerical  error  versus  resolving  grid 
cell  size,  <$jc,  indicating  the  convergence  property  of  the  quantum 
lattice-gas  algorithm  (12)  and  (13)  for  the  Schrodingcr  equation. 
The  data  (black  circles)  are  taken  from  numerical  simulations  with 
grid  sizes  from  L  =  St  up  to  8192f  after  a  single  time  step  T  —  r. 
The  solid  curves  are  best-fit  linear  regression  with  a  slope  of  3.48 
and  5.45  for  the  models  defined  by  (12)  and  (13),  respectively. 
These  results  demonstrate  third-order  and  fourth-order  convergence 
in  space  for  the  two  models,  respectively. 
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3.  Adding  an  external  scalar  potential 


It  is  possible  to  model  an  external  potential  by  ap¬ 
plying  a  local  phase  change  to  the  system  wave  func¬ 
tion  [13,14] 

t)  -»  e”,V(jf)5/i A(jc,  t).  (31) 

The  effect  of  this  phase  change  is  to  alter  the  finite 
difference  equation  (20)  as  follows 

+  e-iV{x'-')S,xHXl-i,t0)] 

-^-[e-iV{XM)S,iHxi+2,to) 

+  e-'V{xi-2)S,rHxi-2Jo)].  (32) 

If  we  expand  the  potential  terms  in  the  arguments  of 
the  exponentials 


V(xt+\)8t  =  V(xi)8t  +  8t8x 


dV(x) 


dr 


X—Xl 


+  0(8t8x2) 
(33) 


we  see  that  we  can  neglect  the  second  term  on  the 
RHS  because  of  diffusive  ordering  8t8x  ~  8; c3  since 
we  need  to  keep  terms  only  to  order  8x2.  Therefore,  in 
the  continuum  limit  (32)  is  well  approximated  by 


8) 

=  -^e-iV{x,)S,nx,,to) 

+  c-'V(x,)S,[t(x,+uto)  +  iHx,-Ut0)] 

-  ^-e-'V(xi)S‘[t(xl+2,  t0)  +  t(xi-2,  fo)]- 

(34) 

Now  multiplying  through  by  elV,(j/)$'  and  expanding 
the  LHS  to  order  St 2  we  have  the  following  finite- 
difference  equation: 


[l+iV(jr, )«/]*(*,,  *8) 

1  +  •  , ,  x 

=  -j-^(*Mo) 

+  'l'(xi+i,to)  +  \K*/-Mo) 


1  -i 


[lK*/+2.  to)  +  V^U/-2.  to)]. 


(35) 


In  the  continuum  limit,  this  finite-difference  equation 
represents  the  Schrodinger  wave  equation  with  an 
external  potential  term 


9t/c(jc,  t) 

Ft 


+  (D(8t2) 


i  Sx2  d2\fr(x,  t) 

2~S7 


-iV(x)\(r(x,t)  +  0(8x4). 


(36) 


To  confirm  the  validity  of  (36)  we  perform  the  follow¬ 
ing  numerical  simulations  that  yield  results  that  can  be 
checked  against  analytical  predictions: 


(1)  Harmonic  oscillation  of  a  displaced  Gaussian 
wave  packet  in  a  parabolic  potential. 

(2)  Quantum  tunneling  through  a  potential  barrier. 


3. 1.  Harmonic  oscillator 


The  first  numerical  test  presented  here  is  the  simu¬ 
lation  of  the  behavior  of  a  wave  packet  in  an  external 
parabolic  potential.  This  is  the  well-known  problem 
of  the  linear  harmonic  oscillator.  Schrodinger  analyti¬ 
cally  calculated  the  exact  time-dependent  solution  for 
the  evolution  of  a  Gaussian  packet  that  is  displaced 
by  a  distance  a  from  its  central  ground  state  in  a  par¬ 
abolic  potential  well  of  the  form  V(jc)  =  \Kx2.  The 
initial  wave  function  is 

/y  1  /2  ,  , 

^U,0)  =  ^T7Ie-"2("-fl)2/2,  (37) 

71 XH 

where  a  =  {mK /h2)x^  is  the  width  of  the  packet 
and  coc  =  ( K/m )1^2  is  the  angular  frequency  of  the 
classical  harmonic  oscillator  [21].  The  exact  time- 
dependent  solution  for  the  probability  density  is  the 
following: 


\n*.t)  |2 


a  -a2  (x -a  cos  (Oct)2 /2 

7T ^ 


(38) 


A  derivation  of  the  result  (38)  is  also  presented  by 
Schiflf  [22]. 

To  test  the  quantum  lattice  gas  algorithm  against 
(38)  we  used  a  periodic  lattice  with  L  =  256 1  nodes. 
The  initial  Gaussian  packet  is  displaced  to  the  right 
of  the  center  of  the  grid  by  32  lattice  nodes  and  so 
is  initially  located  at  jto  =  160f  as  shown  in  Fig.  3. 
With  h  =  1  and  m  =  1,  the  classical  time  period  is 
Tc  =  2tt /cdc  =  1987r.  So  letting  the  simulation  run 
for  1000  iterations  allows  the  packet  to  the  other  side 
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of  the  potential  well  near  position  x  =  96 1  as  demon* 
strated  in  Fig.  3. 

The  simulation  was  run  for  a  total  of  6000  time 
steps  and  the  location  of  the  peak  of  the  Gaussian  wave 
packet  was  recorded  every  lOOr  time  steps.  This  data 
is  plotted  in  Fig.  4.  The  location  of  the  peak  oscillates 
in  time  as  expected.  Overplotted  on  this  numerical  data 
is  the  exact  solution  for  the  oscillation  a  cos  coct  +  jcq 


Fig.  3.  Time  evolution  of  a  Gaussian  packet  initially  displaced  by 
a  =  32 1  lattice  sites  from  the  center  of  a  parabolic  potential  well 
with  K  =  10-5.  The  width  of  the  packet  is  a  =  14.4f.  The  time 
development  of  the  Gaussian  packets  over  plotted  in  succession 
where  the  jc-axis  is  the  position  on  a  L  —  256 1  node  lattice  and 
the  y-axis  is  the  probability  density  \\J/(x,t)\^.  The  red  curve  is 
the  parabolic  potential.  The  h  =  1  and  m  —  1,  the  time  period  of 
the  oscillation  is  Tc  =  ~  =  1986.92r.  A  total  of  ten  profiles  are 

over  plotted  corresponding  to  time  /  =  0,  100r,200r _ _  lOOOr, 

which  is  approximately  half  of  the  oscillation  time  period,  so  the 
packet  is  seen  to  “swing”  to  the  other  side  of  the  potential  well  while 
maintaining  a  fixed  shape  as  analytically  predicted. 


Time  Step  Iteration  ( X ) 

Fig.  4.  A  comparison  between  the  analytical  and  numerical  predic¬ 
tions  of  the  location  of  an  oscillating  Gaussian  packet  in  a  harmonic 
parabolic  potential  well.  The  solid  curve  is  the  analytical  prediction 
and  the  black  circles  are  the  numerical  data  taken  from  the  quan¬ 
tum  lattice  gas  simulation  presented  in  Fig.  3.  In  the  simulation, 
the  packet  is  initially  displaced  32  lattice  units  from  the  center  of 
the  grid  at  lattice  node  128  for  a  periodic  system  with  a  total  of 
L  =  256 1  nodes.  The  numerical  predictions  are  in  excellent  agree¬ 
ment  with  the  exact  analytical  solution. 


and  the  agreement  between  the  analytical  solution  and 
the  numerical  data  is  excellent. 

3.2.  Scattering  off  a  potential  barrier 

The  next  numerical  test  of  the  quantum  lattice  gas  is 
to  simulate  the  well-known  case  of  quantum  tunneling 
through  a  constant  potential  barrier  of  width  a.  That 
is,  V(x)  =  Vo  for  0  ^  x  ^  a  and  V(jt)  =  0  otherwise. 
The  initial  wave  function  is  a  Gaussian  packet  with  net 
momentum  to  the  right 

=  (39» 

where  p  is  the  momentum  parameter.  We  choose  the 
mean  kinetic  energy  of  the  packet  to  be  equal  to  the 
constant  energy  level  of  the  potential  barrier  \  p2  = 
Vo*  In  this  case,  the  packet  tunnels  through  the  barrier 
but  the  sum  of  the  transmission  and  reflection  proba¬ 
bilities  are  less  than  one  because  there  is  a  resonance 


Fig.  5.  A  sequence  of  snap  shots  of  the  time  evolution  of  a 
packet  that  is  incident  from  the  left  on  to  a  potential  barrier 
where  the  mean  kinetic  energy  of  the  packet  equals  the  energy 
of  the  barrier.  The  x-axis  is  the  lattice  position  and  the  y-axis  is 
the  probability  density.  The  iteration  time  step  for  each  frame  of 
the  sequence  is  labeled  in  the  upper  left  comers.  The  simulation 
was  run  on  a  periodic  grid  of  size  L  =  4000f  for  a  total  of 
20,000  time  steps.  The  width  of  the  incident  packet  was  set  to 
cr  =  0.035L  =  140f  and  the  initial  momentum  parameter  was  set 
to  p  =  0. 1  in  units  where  t=  1 ,  r  =  1  and  m  =  1 .  The  width  of  the 
barrier  was  set  to  a  =  0.064/.  =  256f .  As  expected  the  numerical 
simulation  clearly  demonstrates  the  resonance  effect  where  there  is 
a  non-zero  probability  of  the  particle  to  be  trapped  within  the  barrier 
itself. 
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effect  where  the  particle  is  also  trapped  inside  of  the 
barrier.  This  effect  is  observed  in  the  numerical  simu¬ 
lation  shown  in  Fig.  5. 


4.  Two  fermionic  particles 


"maxr  /  x  \ 
(p(xj)  =  +  an  cos^2h7t  —  J 

4-  bn  s'\n^2nn 


iEnt  /  is 


(44) 


The  efficiency  of  the  quantum  algorithm  (12)  be¬ 
comes  evident  when  it  is  used  to  simulate  the  dynam¬ 
ics  of  multiple  quantum  particles.  The  case  of  multiple 
quantum  particles  is  still  handled  by  the  same  evolu¬ 
tion  operator  E  that  we  tested  for  the  single  particle 
case.  The  particular  sequence  and  number  of  quan¬ 
tum  gate  operations  remains  fixed,  independent  of  the 
number  of  particles  to  be  simulated.  The  only  differ¬ 
ence  is  how  the  system  wave  function  is  initialized. 

In  this  section,  for  pedagogical  reasons,  we  will 
consider  the  case  of  simulating  two  free  quantum  par¬ 
ticles.  The  approach  we  use  in  this  case  can  be  directly 
generalized  to  the  many-particle  case. 

To  begin  with  we  write  the  Schrodinger  wave 
equation  for  two  free  quantum  particles 

-d}/f(x,y,t) 

l  h - 

dt 

_  /r_  3V(jc,>’,0  _  P_  d2i/f(x,y,t) 

2m  dx2  2m  dy2 

where  x  and  y  are  the  spatial  coordinates  of  the 
first  and  second  particle,  respectively.  Since  the  wave 
function  is  spatially  separable  as  \J/(x,y,  t)  =  tp{x,  t)  • 
<p(y ,  t)f  the  analytical  solution  to  (40)  is  obtained  in  a 
similar  fashion  to  the  one-body  case  by  computing  the 
Fourier  components  of  the  energy  basis  functions 


ao 


dn 


b 


n 


<p(x,  0)  djc, 


c,  0)cos 


^2/z7T  cLc, 


<p(jc,  0)  sin (  2nn—  )  dv 


r> 


(41) 


(42) 


(43) 


The  energy  eigenvalues  are  still  given  by  (27)  and  the 
time-dependent  single-particle  solution  is 


which  is  basically  the  same  as  (28)  except  that  we 
had  to  add  the  sin  term  because  with  two  particles  the 
wavefunction  is  not  even,  as  is  (24),  for  example.  We 
shall  test  the  time  evolution  of  two  Gaussian  packets. 
The  initial  wave  function  in  our  test  is  the  odd  function 


^AexactC* »  y i  t) 


1 

7 1 


\<Pa\ ,<T| 


(x,t)<Pa2,a2(y,t) 


-<Pa\.oi(y,t)<Pai.oi(x,t)\, 


(45) 


where 

gW(*,0  )  =  gi/217r,/4e— U~c)2/(2<r2).  (46) 

The  subscripts  on  the  function  denote  its  de¬ 
pendence  on  the  position  and  width  of  the  individual 
Gaussian  packet.  This  functional  dependence  is  actu¬ 
ally  contained  within  the  form  of  the  coefficients  ao , 
an ,  and  bn  that  depend  on  the  position  and  width  of  the 
Gaussian  packet  in  accordance  with  (41)  through  (43). 
Note  that  given  the  form  of  (45),  inexact (*>*>  0  =  0 
which  satisfies  the  Pauli  exclusion  principle. 


4.1.  Numerical  confirmation 


To  numerically  simulate  the  evolution  of  the  two- 
particle  wave  function  governed  by  (40)  using  quan¬ 
tum  algorithm  (12)  we  must  use  a  new  computational 
formulation  to  implement  our  algorithm.  The  finite- 
difference  equation  implementation  that  we  used  in 
Section  2.4,  in  the  single-particle  case,  cannot  be  di¬ 
rectly  applied  in  the  two-particle  case  to  each  particle 
individually  because  it  does  not  allow  for  the  possibil¬ 
ity  when  the  particles  are  quantum  mechanically  en¬ 
tangled.  In  general,  this  will  be  the  case  when  there 
is  an  interaction  between  the  particles.  Therefore,  we 
shall  use  an  implementation  that  can  handle  the  most 
general  situations  involving  correlated  particles  and 
one  that  naturally  scales  to  handle  an  arbitrarily  large 
number  of  particles  in  the  system. 
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We  shall  represent  the  basic  quantum  gate  opera¬ 
tions  in  terms  of  the  fermionic  creation  and  annihila¬ 
tion  operators  in  the  number  representation,  denoted 
al  and  aa ,  respectively,  and  use  this  approach  as  the 
basis  for  a  general  computational  formulation  applica¬ 
ble,  in  particular,  to  our  algorithm  and,  in  general, 
to  any  quantum  algorithm.  Acting  on  a  system  of  Q 
qubits,  al  and  aa  create  and  destroy  a  particle  occu¬ 
pancy  encoded  in  the  ath  qubit 


al\n\ ... na...nQ ) 


aa|«l  ...na...nQ) 


0  na  —  1, 

|«l  ...1  ...«e)  na  =  0, 

(47) 

\n\...0...nQ)  na  =  1, 

0  na  —  0. 

(48) 


The  fermionic  creation  and  annihilation  operators 
satisfy  the  anti-commutation  relations 

{oa,  ajj]  =  8afl< 

(aa,ap}  =  0,  (49) 

( ala\)  =  0. 

The  number  operator  na  =  alcia  has  eigenvalues  of 
1  or  0  in  the  number  representation  when  acting  on 
a  pure  state,  corresponding  to  the  orth  qubit  being  in 
state  |1)  or  |0),  respectively. 

The  square-root-of-swap  gate  (5)  acting  on  the  on¬ 
site  qubits  indexed  by  a  and  a  -f  1  can  be  expressed  in 
terms  of  the  creation  and  annihilation  operators  as 

^a,cr  +  l  =  ^  —  Ba+\)  ActyCta+i  —  ^^+1^ 

+  A*(l  -na)na+\  +  1  -ha  -  na+ 1, 

(50) 

where  A  =  5  +  j.  Also,  the  swap  operator  (9)  acting 
between  the  first  qubits  indexed  by  a  and  at 
neighboring  nodes  can  be  expressed  in  terms  of  the 
creation  and  annihilation  operators  as 

Xa.p  =  1  -alap  -c/paa  -ha  -rip.  (51) 

The  quantum  gates  (50)  and  (51)  are  used  to  imple¬ 
ment  the  quantum  lattice  gas  collision  and  streaming 
operations,  respectively  [23]. 

The  basis  state  in  the  two-particle  sector  can  be 
labeled  with  the  binary  encoding  formula  |2Q'-1  + 
2^-1)  where  the  integers  a  and  ft  are  in  the  ranges 


1  ^  a  ^  Q  and  a  -f  1  ^  ^  Q.  The  number  of  basis 

states  in  this  case  is  the  binomial  coefficient  (?).  The 
system  ket  can  then  be  expressed  in  the  two-particle 
sector  as 

Q  Q 

m  =  E  *«./> |2“-,+2*-1>.  (52) 

ct=\  fi=a+\ 

Since  there  are  two  qubits  per  site,  we  initialize  the 
wave  function  using  (45)  as  follows: 


=  V'exactl 


a  +  1 
2 

l+l 

2 


L  -T  1 
2 


(53) 


where  the  notation  means  the  floor  of  x  and 
where  Q  =  2 L.  The  floor  operation  is  used  so  that 
the  initial  value  of  the  wave  function  at  each  node  is 
divided  evenly  between  each  pair  of  on-site  qubits. 
This  is  on  account  of  definition  (4)  that  allowed  us 
to  have  interfering  possibilities  for  a  single  particle 
to  occupy  a  single  position  on  the  lattice.  Moreover 
in  the  two  particle  case,  still  only  a  single  particle 
can  occupy  a  single  position  because  of  the  form  of 
the  wave  function  (45)  which  is  consistent  with  the 
anti-commutator  relations  (49).  However,  particle  one 
can  interfere  on-site  with  itself  or  with  particle  two, 
or  vice  versa  since  the  particles  arc  indistinguish¬ 
able. 

At  this  point  we  have  described  how  we  implement 
the  two  quantum  gates  used  in  our  algorithm,  how 
we  enumerate  the  basis  states,  and  how  we  initialize 
the  two-body  wave  function  in  this  basis.  The  only 
remaining  issue  left  to  describe  is  how  we  project 
the  two-coordinate  wave  function  \ //(x,y,t)  on  to 
a  single-coordinate  wave  function  \l/(x,t)  that  can 
be  plotted  on  a  single  physical  axis.  Because  of  the 
underlying  lattice  in  our  system,  this  is  straightforward 
to  do  by  summing  out  one  of  the  coordinates  as 
follows: 


L- 1 

'I'ixi,  tn)  =  ^2  iM*/.  ym,  tn)-  (54) 

m=  0 


If  \l/(xiyym,tn)  is  normalized  then  so  is  \f/{x /,/„) 
according  to  (54).  A  comparison  of  the  time  evolution 
of  the  analytical  solution  (45)  and  the  numerical 
solution  (54)  for  a  lattice  with  30  nodes  is  shown  in 


J.  Yepez,  B.  Boghosian  /  Computer  Physics  Communications  146  (2002)  280-294 


292 


5  10  15  20  25  30 

Position 


Fig.  6.  Time  evolution  of  two  fermionic  particles  initialized  as 
Gaussian  packets  overplotted  in  succession  where  the  jt-axis  is  the 
position  on  a  30-node  lattice  in  units  of  the  lattice  spacing  t  and 
the  v-axis  is  the  probability  density  \\J/(xy ,  x-},  t)\2  projected  onto 
the  x |  -axis.  The  solid  curves  are  the  exact  analytical  solution  and 
the  circles  are  the  data  from  the  quantum  lattice-gas  simulation  (the 
initial  wave  function  was  normalized,  therefore  the  area  under  each 
curve  is  one).  The  initial  Gaussian  packets  of  width  cr  =  3f  at  t  =  0 
of  the  first  and  second  particle  is  centered  at  x  =  10 i  and  x  =  20£, 
respectively.  The  dispersion  of  both  packets  is  evident  by  observing 
the  wave  function  at  the  subsequent  times  /  =  7r,  21r,  28r,  35r 
and  42r.  Periodic  boundary  conditions  were  used  and  nmax  =  40 
energy  eigenmodes  were  used  to  generate  the  exact  solutions  at  four 
times  the  resolution  of  the  numerical  solution.  No  time  scale  factor 
was  used  and  there  is  good  agreement  between  the  analytical  and 
numerical  predictions  at  all  later  times  of  the  numerical  simulation 
as  demonstated  by  the  graphs. 

Fig.  6.  Even  with  this  small  lattice,  throughout  the  time 
evolution  of  the  model  run  the  numerical  predictions 
are  in  good  agreement  with  the  predictions  of  the  exact 
solution. 


5.  Conclusion 

We  have  presented  a  quantum  algorithm  that  is 
an  efficient  and  accurate  scheme  for  simulating  the 
time-dependent  evolution  of  a  system  of  quantum 
particles  governed  by  the  non-relativistic  Schrodinger 
wave  equation.  The  scheme  uses  a  quantum  lattice 
gas  system  of  particles  colliding  and  hopping  on  a 
lattice.  The  algorithm  is  efficient  in  the  sense  that  the 
computational  effort  needed  to  simulate  an  arbitrarily 
large  number  of  particles  (within  the  constraint  of 
the  grid  resolution)  exactly  equals  the  computational 
work  needed  to  simulate  a  single  particle,  given  that 
the  algorithm  is  executed  on  a  quantum  computer 
that  remains  phase-coherent  throughout  the  entire 


coarse  of  the  simulation.  However,  a  limitation  does 
exist  on  state  preparation  (i.e.  initialize  the  quantum 
computer’s  memory)  and  we  have  not  argued  here 
that  it  is  possible  to  initialize  the  many-body  wave 
function  in  an  efficient  way.  For  that  matter,  have  we 
also  have  not  argued  that  it  is  possible  to  measure  the 
final  state  (reading  the  quantum  computer’s  memory) 
of  the  computed  wave  function  in  an  efficient  way. 
Nevertheless,  the  quantum  algorithm  presented  here, 
which  is  a  way  of  representing  a  discretized  Feynman 
path  integral,  has  the  useful  feature  that  it  is  explicit  in 
time  where  the  value  of  the  wave  function  at  location  jc 
and  time  /  +  :  depends  only  on  the  previous  values  of 
the  wave  function  at  time  t  in  the  immediate  vicinity 
of  x.  Since  the  algorithm  is  unitary  and  fourth  order 
accurate  in  space,  it  is  useful  even  for  implementation 
on  a  classical  computer. 

We  have  carried  out  a  variety  of  numerical  tests 
proving  that  the  quantum  algorithm  indeed  allows  us 
to  faithfully  reproduce  the  correct  dynamical  behav¬ 
ior  of  a  continuous  and  differentiable  wave  function 
in  the  presence  of  an  external  potential.  However,  the 
total  probability  for  finding  the  quantum  particle  in 
the  system  is  not  exactly  conserved  in  this  quantum 
lattice  gas  model.  One  must  approach  the  continuum 
limit  to  achieve  a  high  degree  of  probability  conser¬ 
vation.  Since  we  have  demonstrated  that  the  quantum 
lattice  gas  model  is  fourth-order  convergent  in  space, 
it  is  always  possible  to  choose  a  grid  resolution  that 
achieves  the  necessary  fixed  numerical  accuracy  re¬ 
quired  by  any  application. 

We  have  also  described  and  carried  out  the  numer¬ 
ical  simulation  of  two  fermionic  particles,  which  are 
non-interacting  except  for  a  quantum  mechanical  ex¬ 
change  force  arising  from  the  anti-commutation  re¬ 
lations.  The  numerical  formalism  used  to  implement 
the  quantum  lattice  gas  algorithm  represents  the  ba¬ 
sic  quantum  gates  in  terms  of  quadratic  products  of 
creation  and  annihilation  operators  in  a  second  quan¬ 
tized  representation.  This  formalism  straightforwardly 
handles  an  arbitrarily  large  number  of  particles.  The 
simulation  is  carried  out  in  the  many-body  sector  (i.e. 
a  Fock  space  where  the  number  of  particles  is  fixed) 
where  all  the  basis  states  are  enumerated  by  a  simple 
binary  encoding  formula.  In  general,  the  wave  func¬ 
tion  is  initialized  using  a  Slater  determinant  so  that 
the  Pauli  exclusion  principle  is  satisfied  and  the  wave 
function  is  odd.  In  all  the  test  cases  (single  free  par- 
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tide,  single  particle  in  an  external  potential,  and  two 
free  fermions)  the  numerical  predictions  agreed  ex¬ 
tremely  well  the  analytical  predictions  and  exact  so¬ 
lutions. 

It  is  important  to  realize  that  it  is  possible  to  push 
the  quantum  lattice  gas  model  into  regimes  where  the 
numerically  predicted  results  are  absolutely  wrong. 
This  occurs  when  the  local  configurations  are  far  from 
local  equilibrium.  Remember  that  local  equilibrium 
exists  on  a  lattice  node  when  the  qubits  at  that  node 
have  identical  phase.  Therefore,  large  gradients  in  the 
macroscopic  profile  of  the  modeled  wave  function 
may  cause  a  large  phase  difference  between  the  on-site 
qubits  eventually  resulting  in  anomalous  large-scale 
behavior  in  the  model.  A  good  example  of  this  occurs 
when  the  momentum  of  a  moving  wave  packet  is  too 
high.  In  this  case,  since  the  real  and  imaginary  parts 
of  the  wave  function  are  sinusoidal,  if  the  momentum 
is  so  high  that  the  wavelength  of  the  traveling  wave 
is  on  the  order  of  the  lattice  cell  size  (A  =  h  / p  ~  f), 
then  after  a  few  times  step  iterations  of  the  algorithm, 
large  phase  differences  in  the  on-site  qubits  occur. 
In  this  case,  the  local  on-site  configuration  will  not 
relax  toward  the  correct  local  equilibrium.  Fortunately, 
the  norm  of  the  modeled  wave  function  will  deviate 
from  unity  in  these  types  of  cases.  Therefore,  it  is 
straightforward  to  test  if  the  model  predictions  are 
non-physical  by  periodically  checking  the  norm  of  the 
wave  function. 

Interaction  potentials  between  particles  can  also  be 
modeled  using  this  quantum  lattice  gas  method  [14]. 
We  still  need  to  perform  simulations  of  a  many-body 
system  with  an  interaction  potential.  Also,  numerical 
tests  in  two  and  three  dimensions  should  be  conducted. 
In  this  case,  it  would  be  straightforward  to  test  the  ad¬ 
dition  of  an  external  vector  potential  using  the  results 
analytically  predicted  by  Policy  [18].  The  simulation 
of  the  dynamical  behavior  of  the  positronium  atom 
would  be  a  reasonable  next  step  for  the  application  of 
the  quantum  lattice  gas  method  to  quantum  mechani¬ 
cal  computational  physics. 
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Appendix  A 


The  full  finite-difference  equation  for  the  quan¬ 
tum  lattice-gas  model  presented  in  this  paper  is  quite 
long.  To  simplify  this  expression,  we  introduce  a  lo¬ 
cal  neighborhood  vector  with  the  following  18  com¬ 
ponents 

rj(xi,t) 

=  (<po(*/,0,<Pi(*/,0, 

W)(*/+i.  0,  <P\  U/+ 1,0,  <po(*/-i>  0,  v  i  (*/- 1,  Of 
W>(*/+2»  0.  <P\ (*/+2.  0.  <Po(*l-2f  Of  <P\  (*/- 2*  Of 
W>(*/+3.  Of  <P\  (*/+3»  Of  <Po(xi-3f  Of  <P\  (xi- 3,  Of 

<poU/+4,  Of  <p\  (*/+4.  Of  <PoU/-4,  Of  <P\  (*/— 4,  0). 

(A.  1 ) 

We  define  the  following  two  coefficient  vectors 


a  =  (—3,  — 3i,  3,  i,  —3,  — 5i,  1,  i,  3,  3i, 

-  1,  i,  1, 3i,  0, 0,  —1,  — i), 

P  =  ( — 3i,  -3,  — 5i,  -3,  i,  3,  3i,  3,  i,  1, 3i, 
1,  i,  — 1,  — i,  —1, 0, 0). 


(A. 2) 


The  microscopic  evolution  equation  (13),  explicitly 
written  out,  has  the  following  protocol  of  operations 


|  *(fie))  =  (S/CS2CS/CS2C) 

X  (^CSlC^CSiCjl^dtt)).  (A. 3) 

The  corresponding  full  finite-difference  equation  can 
be  specified  by  the  following  dot  product  of  these 
vectors 


<Po(*Mi6)  = 


<Pi(xi,t\6)  = 


a  •  q(*;.  to) 
16 

P-  q(*Mo) 
16 


(A. 4) 


(A. 5) 


Note  that  if  to  is  the  initial  time,  then  the  interval 
r  =  /] 6  is  defined  the  update  time  step.  The  finite- 
difference  equation  for  \J/  =  <po  +  <p\  is 


'friXht  16)  = 


(a +  /?)•>?(•*;,  fo) 

16 


(A. 6) 


and  it  has  a  high  degree  of  numerical  accuracy  as 
indicated  in  Fig.  2. 
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Abstract 

I  review  the  differences  between  classical  and  quantum  systems,  emphasizing  the  connection  between  no-hidden  variable 
theorems  and  superior  computational  power  of  quantum  computers.  Using  quantum  lattice  gas  automata  as  examples,  I  describe 
possibilities  for  efficient  simulation  of  quantum  and  classical  systems  with  a  quantum  computer.  I  conclude  with  a  list  of  research 
directions.  ©  2002  Published  by  Elsevier  Science  B.V. 

PACS:  03. 67. Lx;  05.10.-a 
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1.  Introduction 

There  are  two  paths  towards  quantum  computing: 
one  is  teleological  and  the  other  is  practical.  The  tele¬ 
ological  path — described  35  years  ago  in  the  proph¬ 
esy  known  as  Moore’s  Law  [1] — leads  down  through 
smaller  and  smaller  device  sizes  where  quantum  ef¬ 
fects  become  wilder  and  wilder.  Eventually,  rather 
than  domesticating  them  for  classical  computation,  ex¬ 
perimental  physicists  and  engineers  believe  they  will 
be  able  to  preserve  them  for  quantum  computation. 
The  practical  path,  on  the  other  hand,  is  paved  with 
the  desire  to  solve  specific  problems  efficiently.  In  an 
amusing  role-reversal,  it  is  theoretical  physicists,  com¬ 
puter  scientists,  and  mathematicians  who  follow  this 
path.  The  first  steps  along  it  were  taken  20  years  ago 
by  Feynman,  who  suggested  that  since  quantum  sys¬ 
tems  seem  to  be  very  hard  to  simulate  on  a  classical 
computer,  perhaps  they  could  be  simulated  more  ef- 
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ficiently  on  a  quantum  computer  [2].  More  recently, 
Deutsch  [3],  Jozsa  [4],  Simon  [5],  Shor  [6],  Grover  [7] 
and  others  have  noted  that  a  quantum  computer  could 
solve  classical  problems  as  well.  In  this  primarily  ped¬ 
agogical  paper  1  describe  some  of  the  steps  which 
have  been  taken  along  this  practical  path,  and  spec¬ 
ulate  about  some  steps  further  along  it. 

2.  Simulating  quantum  systems  classically 

Let  me  begin  by  reviewing  the  reasons  quantum 
systems  are  believed  to  be  hard  to  simulate  on  clas¬ 
sical  computers.  Traditionally  these  are  known  as  ‘no 
hidden  variable  theorems’.  Each  is  a  statement  that 
no  classical  model  with  specified  constraints  can  re¬ 
produce  quantum  mechanical  results.  Consideration 
of  two  of  them,  the  Gleason/Kochen-Specker  theo¬ 
rem  [8,9]  and  Bell’s  theorem  [10],  reveals  both  their 
heuristic  power  and  their  weaknesses. 

In  1957  Gleason  proved  that  for  Hilbert  spaces  of 
dimension  at  least  3,  any  non-negative  measure  on 
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states  which  is  quantum  mechanical  (i.e.  for  any  or¬ 
thogonal  basis  {£/}  the  measure  sums  to  1)  must  de¬ 
rive  from  a  density  matrix  [8].  In  1967  Kochen  and 
Specker  made  the  contradiction  with  a  classical  hidden 
variable  model  more  explicit  [9]:  They  constructed  a 
finite  set  of  unit  vectors  in  IR3  with  the  property  that 
every  attempt  to  assign  values  0  or  1  to  each  vec¬ 
tor  satisfying  the  condition  that  in  each  orthogonal 
triple  two  vectors  get  1  and  the  third  gets  0  must  fail. 
That  is,  no  classical  ‘hidden  variable’  can  be  assigned 
to  pre-determine  which  outcome  of  each  of  some  fi¬ 
nite  set  of  complete  measurements  of  the  spin-squared 
of  a  spin-1  particle  will  be  observed  (since  the  spin 
eigenvalues  are  { — 1,0,  4- 1 } ,  two  of  the  spin-squareds 
are  1  and  one  is  0  for  any  complete  measurement). 
Such  a  hidden  variable  would  be  non-contextual ,  in 
the  sense  that  its  value  on  each  vector  would  specify 
the  spin-squared  observed  for  that  measurement,  inde¬ 
pendently  of  which  complete  measurement  including 
it  is  performed.  One  can  argue,  however,  that  noncon- 
textuality  is  too  strict  a  condition  to  place  on  hidden 
variables — perhaps  the  results  of  measurement  should 
depend  on  hidden  variables  inherent  in  the  measur¬ 
ing  device,  which  might  differ  for  each  complete  mea¬ 
surement  [11].  Furthermore,  the  measurements  must 
be  exactly  along  the  vectors  Kochen  and  Specker  con¬ 
structed,  but  from  a  computational  complexity  per¬ 
spective,  infinite  precision  is  suspect  [12] — and  in  fact 
one  can  show  that  without  additional  assumptions  one 
cannot  prove  a  Kochen-Specker  theorem  using  only 
finite  precision  measurements  [13]. 

Both  of  these  weaknesses  are  absent  from  Bell’s 
theorem  [10].  He  proved  that  the  results  of  local 
measurements  on  specific  states  of  pairs  of  spin-^ 
particles,  i.e.  vectors  in  C2  <8>C2,  cannot  be  reproduced 
by  any  local,  classical  hidden  variable  model.  Here 
‘local’  means  restricted  to  individual  particles.  This 
result  is  robust  against  measurement  imprecision,  and 
locality  of  the  hidden  variables  seems  justified  on 
physical  grounds — the  finite  speed  of  light  and  the 
locality  of  physical  interactions.  In  fact,  these  are  the 
same  grounds  upon  which  we  base  our  models  of 
computation:  At  each  timestep  a  classical  or  quantum 
Turing  machine  changes  only  the  state  of  the  head  and 
the  symbol  written  on  the  tape  cell  where  the  head  is 
located  [3,14];  it  does  not  make  non-local  changes  of 
all  the  cells  of  the  tape  simultaneously,  for  example. 


The  states  for  which  Bell’s  theorem  rules  out  classi¬ 
cal  hidden  variables  are  entangled ,  i.e.  ones  for  which 
the  state  of  multiple  particles  cannot  be  described  as 
the  product  of  states  for  each  particle  individually. 
Since  this  is  true  for  all  but  a  set  of  measure  0  in  the 
space  of  all  pure  states,  Bell’s  theorem  and  its  gener¬ 
alizations  (see  [15]  for  a  recent  survey)  indicate  that 
most  quantum  states  cannot  even  be  described  by  rea¬ 
sonable  (in  the  sense  of  local)  classical  models.  This 
is  a  more  subtle  problem  than  simply  the  large  size  of 
the  state  space,  which  we  consider  next. 

3.  Dynamics 

The  dimension  of  the  Hilbert  space  describing  the 
state  of  a  system  of  multiple  particles  grows  exponen¬ 
tially  in  the  number  of  particles:  2n  for  n  spin-  \  par¬ 
ticles,  for  example.  This  exponential  explosion,  how¬ 
ever,  is  not  enough  to  preclude  classical  simulation. 
Consider  a  classical,  probabilistic  lattice  gas.  On  a  ho¬ 
mogeneous  one-dimensional  lattice  of  size  n  there  are 
4n  basis  states  s, ,  since  each  lattice  site  can  be  occu¬ 
pied  by  no  more  than  1  particle  with  each  of  the  two 
possible  velocities.  A  general  state  s  is  a  convex  com¬ 
bination: 

4"  — l  4n  — l 

S=J2  PiSi  with  X  Pi  =  1  ’  Pi  ^  °- 

i=0  /=0 

Evolving  the  whole  state,  i.e.  the  probability  distrib¬ 
ution,  is  therefore  an  exponentially  difficult  problem 
in  the  size  of  the  lattice.  Nevertheless,  such  lattice 
gas  models  are  used  regularly  (see,  e.g.,  [16]).  But 
one  does  not  evolve  the  whole  probability  distribution. 
Rather,  one  samples  it,  by  evolving  a  single  s/  to  a  sin¬ 
gle  si  at  the  next  timestep,  using  some  random  num¬ 
ber  generator.  Multiple  runs  sample  the  final  distribu¬ 
tion.  A  quantum  lattice  gas  automaton  (QLGA,  which 
I  will  describe  in  more  detail  in  §4)  is  also  described  at 
each  timestep  by  a  vector  in  a  space  with  basis  — 
where  |*>  is  the  standard  Dirac  notation  for  vectors  in 
Hilbert  space  [17]: 

4n  —  1  4'1  — 1 

1 1)  =  ^  aj\si)  with  ^|fl(|2  =  l,  a;  6  C. 

i=0  i=0 

Evolving  the  QLGA  state  has  classical  computational 
complexity  comparable  to  evolving  the  whole  state  of 
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the  probabilistic  LGA.  But  in  the  quantum  case,  this 
cannot  be  reduced  by  sampling  individual  histories: 
each  has  a  complex  amplitude  so  the  histories  with 
each  given  final  state  interfere. 

Thus  interference  seems  to  be  the  phenomenon 
which  makes  quantum  dynamics  hard  to  simulate  clas¬ 
sically.  In  fact,  although  the  multi-particle  structure  of 
a  system  is  important,  entanglement  per  se  seems  to 
be  less  relevant:  In  liquid  state  NMR  quantum  com¬ 
puting  experiments  [18],  for  example,  the  state  is  not 
entangled  at  any  timestep  (more  precisely,  since  the 
system  is  in  a  mixed  state — a  convex  combination  of 
pure  states — the  state  is  separable)  [19].  Nevertheless, 
it  seems  to  be  difficult  to  construct  a  reasonable  lo¬ 
cal  hidden  variable  model  for  the  dynamics  [20],  i.e. 
the  dynamics  seems  difficult  to  simulate  classically. 
To  make  this  more  than  heuristic,  however,  we  would 
need  a  dynamical  Gleason/Kochen-Specker/Bell-type 
theorem  which  applies  even  for  evolution  through  a  se¬ 
quence  of  unentangled  states.  Perhaps  some  hint  about 
a  way  to  do  this  may  be  found  in  Laflamme’s  response 
to  the  separability  criticism  of  NMR  quantum  compu¬ 
tation  [21]. 

Of  course,  some  demonstrations  of  the  absence  of 
classical  models  for  quantum  dynamics  already  exist. 
These  are  more  commonly  known  as  quantum  algo¬ 
rithms  for  oracle  problems;  since  each  consists  of  a 
sequence  of  unitary  operations,  they  are  dynamical  re¬ 
sults.  Grover’s  quantum  search  algorithm,  for  exam¬ 
ple,  solves  the  problem  of  identifying  a  e  {0,  1}"  given 
an  oracle  which  responds  to  a  query  x  e  {0,  1 }"  by  re¬ 
turning  8xay  using  only  0(\/2" )  queries  [7].  Classi¬ 
cally,  any  algorithm  would  require  0(2")  queries.  For 
n  >  2,  the  state  is  entangled  at  every  timestep  after 
the  first  [22].  Possibly  more  to  the  point  is  Bernstein 
and  Vazirani’s  algorithm  which  solves  the  problem  of 
identifying  a  e  {0,  1 }"  given  an  oracle  which  responds 
to  a  query  x  e  (0,  1}"  by  returning  x  •  a  mod  2,  us¬ 
ing  only  1  quantum  query  [23].  Classically,  any  al¬ 
gorithm  would  require  O(az)  queries.  And  this  quan¬ 
tum  algorithm  works  without  creating  entanglement  at 
any  timestep  [24].  These  results  suggest  that  while  a 
theorem  on  the  impossibility  of  efficient  classical  sim¬ 
ulation  of  quantum  dynamics  may  exist,  it  will  have 
to  count  all  the  elementary  operations,  not  just  the 
queries,  which  will  presumably  make  it  more  difficult 
to  prove. 


4.  Quantum  simulations 

In  §§1-3  I  have  tried  to  explain  the  heuristic  that 
classical  simulation  of  quantum  systems  is  difficult, 
while  noting  what  remains  to  be  proved  to  make 
such  a  claim  rigorous.  Now  let  us  consider  Feynman’s 
proposed  solution:  simulation  with  quantum  comput¬ 
ers  [2].  The  standard  model  of  quantum  computa¬ 
tion  allows  polynomially  many  local  (i.e.  acting  non- 
trivially  on  only  1  or  2  qubits)  gate  operations  [25]. 
This  is  a  reasonable  model  since  in  principle  it  can 
be  realized  by  a  quantum  system  with  a  local  Hamil¬ 
tonian.  Feynman’s  proposal  has  been  verified  in  this 
model  for  quantum  systems  defined  by  local  Hamil¬ 
tonians  [26-31].  More  exotic  quantum  systems  can 
also  be  simulated  efficiently  with  a  standard  quantum 
computer:  Fractional  quantum  Hall  systems,  for  ex¬ 
ample,  have  Hamiltonians  which  vanish  on  the  phys¬ 
ical  states;  the  only  nontrivial  unitary  transformations 
have  global  (topological)  origin.  Nevertheless,  Freed¬ 
man,  Kitaev  and  Wang  have  shown  that  such  topologi¬ 
cal  quantum  field  theories  can  be  simulated  efficiently 
with  a  standard  quantum  computer  [32]. 

A  particularly  simple  architecture  for  a  quantum 
computer  is  a  QLGA  [33].  Although  I’m  not  aware 
of  a  demonstration  that  classical  LGA  are  capable  of 
universal  computation,  their  similarity  to  the  reversible 
billiard  ball  model  of  Margolus  [34]  suggests  that 
they  may  be;  since  QLGA  specialize  to  deterministic 
LGA,  they  would  be  also.  Whether  they  can  efficiently 
(i.e.  with  polynomial  overhead)  simulate  quantum  gate 
arrays  is,  I  believe,  also  an  open  question.  In  the  other 
direction,  QLGA  can  be  simulated  efficiently  on  a 
standard  quantum  computer,  but  have  theoretical  and 
possibly  practical  advantages:  They  directly  simulate 
quantum  systems  and  are  possibly  more  easily  realized 
experimentally  than  arbitrary  quantum  gate  arrays. 

The  possible  configurations  for  each  particle  on 
a  one-dimensional  lattice  L  are  labeled  by  pairs 
(jt,a)  €  L  x  {±1},  where  x  is  the  position  and  a 
the  velocity.  A  classical  lattice  gas  evolution  rule 
consists  of  an  advection  stage  (x,ct)  (x  +  or,  a), 
followed  by  a  scattering  stage.  Each  particle  in  a 
QLGA  [33]  exists  in  states  which  are  superpositions 
of  the  classical  states:  \x//)  =  £ a),  where 
1  =  (^|t/t)  =  ^2  i'x,a'l/x,ce  The  evolution  rule  must  be 
unitary;  the  most  general  with  the  same  form  as  the 
classical  rule  is: 
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^2'l'x.a\x,a)  Y^'I'x.aU  +a,a) 

scatter  , 

1 — »  /  ,Yx,gS(m'\x  +  ct,a  ), 

where  the  scattering  matrix  is 

S  —  (  cosm  i  sin  am  \ 

\isinm  cosm  J  ’ 

Fig.  1  illustrates  this  quantum  evolution:  at  m  —  0 
it  specializes  to  the  classical  deterministic  lattice  gas 
rule.  The  Ax  =  At  -»  0  limit  of  this  discrete  time 
evolution  is  the  Dirac  equation  [33];  the  Ax2  =  At  — ► 
0  limit  is  the  Schrodinger  equation  [35]. 

This  QLGA  model  can  be  extended  to  include 
multiple  particles  with  a  unitary  two  particle  scattering 
rule 

|a,  a,  x ,  —a)  i->»  e1,y|jt,  a ,  a,  —a) 

shown  in  Fig.  1.  With  these  rules  the  1 -dimensional 
QLGA  discretizes  the  quantum  field  theory  described 
by  the  (1  4-  l)-dimensional  massive  Thirring  model 
[33,36].  These  rules  also  preserve  the  symmetry  (i.e. 
bosonic  or  fermionic)  of  the  wave  function  under  par¬ 
ticle  exchange  [37].  The  QLGA  rules  can  be  general¬ 
ized  to  discretize  the  multi-particle  Schrodinger  equa¬ 
tion  in  arbitrary  dimensions  [35];  it  seems  more  diffi¬ 
cult,  however,  to  create  QLGA  which  discretize  rela¬ 
tivistic  evolution  in  higher  dimensions  [38]. 

The  fact  that  the  QLGA  rules  are  homogeneous, 
i.e.  the  same  at  each  lattice  site  and  at  each  timestep, 
suggests  that  they  might  be  easier  to  implement  than 
general  quantum  gate  arrays  which  are  not.  Possi¬ 
ble  physical  systems  in  which  they  might  be  imple¬ 
mented  include  crystals — as  originally  proposed  by 
Feynman  [2]  and  more  recently  in  the  context  of  solid 
state  NMR  [39] — or  optical  lattices  [40].  A  detailed 
proposal  for  physical  implementation  in  such  systems 
could  motivate  experimental  work  towards  realization 
of  QLGA. 


exp  is 

•— *  •  «— •  - -  •  «  ♦  •  • 

Fig.  1.  The  general  evolution  rules  for  the  one-dimensional  QLGA. 


5.  Simulating  classical  systems 

In  the  previous  sections  we  have  seen  that  there  are 
quantum  algorithms  to  efficiently  simulate  multipar¬ 
ticle  quantum  systems  which  seem  to  be  difficult  to 
simulate  classically.  Since,  as  I  noted  in  the  Introduc¬ 
tion,  there  are  efficient  quantum  algorithms  to  solve 
classical  problems,  a  natural  question  is  whether  quan¬ 
tum  computers  can  simulate  classical  systems  [41,42]. 
Assuming  quantum  mechanics  is  a  correct  description 
of  the  world,  the  existence  of  a  classical  description 
for  macroscopic  physics  means  that  quantum  com¬ 
puters  can  simulate  classical  physics  with  constant 
overhead — although  the  constant  factor  may  be  some¬ 
thing  like  1023,  i.e.  a  number  of  quantum  degrees  of 
freedom  sufficiently  large  that  subsystems  decohere 
and  can  be  identified  as  the  classical  objects  to  be  sim¬ 
ulated. 

Can  we  do  better?  That  is,  could  there  be  quantum 
speedups  for  classical  physics?  Yepez  has  proposed 
that  the  answer  is  ‘yes’.  Using  a  “Type  11“  quantum 
computer  in  which  the  state  is  measured,  locally,  af¬ 
ter  each  timestep  and  then  reset  using  a  lattice  Boltz¬ 
mann  rule  [43].  A  model  like  this  can  achieve  at  most  a 
constant  speedup,  corresponding  to  reduced  computa¬ 
tional  cost  for  local  evolution.  In  practice,  of  course,  a 
large  constant  improvement  can  be  tremendously  use¬ 
ful,  but  perhaps  it  is  possible  to  do  better.  More  pre¬ 
cisely,  using  a  standard  quantum  computer,  can  classi¬ 
cal  systems  be  simulated  more  efficiently  than  is  pos¬ 
sible  classically?  Lidar  and  Biham  have  shown  that 
the  answer  to  this  question  is  also  ‘yes’,  for  the  non- 
dynamical  problem  of  sampling  the  ground  state  dis¬ 
tribution  of  a  spin  glass  [41]. 

There  are  also  QLGA  results  which  suggest  that 
certain  aspects  of  classical  dynamics  can  be  simulated 
more  efficiently  quantum  mechanically.  Consider  clas¬ 
sical  diffusion  of  a  particle  in  a  linear  potential,  as 
shown  in  Fig.  2.  A  discrete  model  for  the  evolution 
is  a  biased  random  walk,  with  prob(AA)  oc  e-VVA*. 
The  results  of  simulating  the  evolution  with  a  classical 
(probabilistic)  LGA  are  shown  in  Fig.  3.  The  average 
position  of  the  particle  satisfies 

(x(t))-[x(0))oc-VVt. 

In  QLGA  the  evolution  rules  can  be  modified  to 
include  a  potential  by  incorporating  an  a  dependent 
phase  multiplication,  i.e.  at  each  timestep  [44], 
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Fig.  2.  The  probability  distribution  for  the  position  of  a  classical 
particle  after  diffusing  in  a  linear  potential.  The  particle  was  initially 
at  the  origin. 

t 


Fig.  3.  The  evolution  of  expected  position  for  the  random  walk 
model  of  a  classical  particle  diffusing  in  a  linear  potential. 

which  one  might  imagine  implementing  in  a  physical 
system  with  an  applied,  spatially  varying  magnetic 
field,  for  example.  Fig.  4  shows  the  result  of  a  QLGA 
simulation  with  a  linear  potential.  Now  the  average 
position  of  the  particle  approximately  satisfies 

(. x(t))-(x(0))<x-VVt 2. 

That  is,  this  quantum  system  simulates  the  evolution 
of  the  average  position  of  a  classical  particle  diffus¬ 
ing  in  a  linear  potential  quadratically  faster  than  does 
the  classical  simulation  shown  in  Fig.  3.  I  must  em¬ 
phasize  that  it  is  only  the  average  position  which  is 
being  simulated  accurately,  not  the  whole  probabil¬ 
ity  distribution.  Furthermore,  the  quadratic  speedup 
only  holds  on  timescales  t  27T/VV.  On  longer 
timescales  the  evolution  is  periodic  [45].  Nevertheless, 
this  very  simple  example  suggests  that  efficient  simu¬ 


Fig.  4.  The  evolution  of  expected  position  for  the  QLGA  model  of  a 
quantum  particle  subject  to  a  linear  potential. 

lation  of  more  complicated  classical  dynamics  may  be 
possible. 

6.  Conclusion 

In  conclusion,  let  me  reiterate  the  open  questions 
discussed  in  this  paper: 

Is  there  a  proof  that  (some)  quantum  dynamics  is 
difficult  to  simulate  classically?  Can  it  be  difficult 
even  when  the  state  is  unentangled  (separable)  at 
each  timestep? 

In  case  QLGA  become  a  practical  architecture  for 
quantum  computers,  can  they  simulate  the  standard 
model  of  quantum  computation  with  no  more  than 
polynomial  overhead? 

What  are  possible  physical  implementations  of 
QLGA? 

What  are  the  correct  QLGA  models  for  relativistic 
quantum  systems  in  more  than  1  spatial  dimension? 

and  most  importantly, 

Are  there  quantum  algorithms  which  speed  up  the 
simulation  of  classical  physics? 

Positive  answers  to  this  last  question  will  broaden  the 
possible  uses  for  a  quantum  computer  and  help  justify 
the  immense  commitment  of  resources  which  seems 
likely  to  be  necessary  to  develop  a  scalable  one. 
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Abstract 

The  real-time  probabilistic  simulation  of  quantum  systems  in  classical  computers  is  known  to  be  limited  by  the  so-called 
dynamical  sign  problem,  a  problem  leading  to  exponential  complexity.  In  1981  Richard  Feynman  raised  some  provocative 
questions  in  connection  to  the  “exact  imitation’’  of  such  systems  using  a  special  device  named  a  “quantum  computer’’.  Feynman 
hesitated  about  the  possibility  of  imitating  fermion  systems  using  such  a  device.  Here  we  address  some  of  his  concerns  and,  in 
particular,  investigate  the  simulation  of  fermionic  systems.  We  show  how  quantum  computers  avoid  the  sign  problem  in  some 
cases  by  reducing  the  complexity  from  exponential  to  polynomial.  Our  demonstration  is  based  upon  the  use  of  isomorphisms 
of  algebras.  We  present  specific  quantum  algorithms  that  illustrate  the  main  points  of  our  algebraic  approach.  ©  2002  Published 
by  Elsevier  Science  B.V. 


1.  Introduction 

Because  of  recent  exciting  algorithms,  like  the  fac¬ 
toring  algorithm  of  Shor  [1]  and  the  search  algorithm 
of  Grover  [2],  that  solve  difficult  problems  on  a  quan¬ 
tum  computer  using  algorithms  that  would  be  im¬ 
practical  on  a  classical  computer,  it  is  easy  to  over¬ 
look  that  the  original  proposals  for  quantum  comput¬ 
ers  were  for  the  purpose  of  solving  quantum  physics 
problems  [3].  People  like  Feynman  [3]  focused  on  the 
extent  to  which  such  a  computer  could  imitate  a  spe¬ 
cific  physical  process,  suggesting  in  part  that  quantum 
problems  were  inherently  too  complex  for  a  classical 
computer  [3]. 

The  obvious  difficulty  with  deterministically  solv¬ 
ing  a  quantum  many-body  problem  (of  fermions  or 
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bosons)  on  a  classical  computer  is  the  exponentially 
large  basis  set  needed  (i.e.  the  dimension  of  its  Hilbert 
space  grows  exponentially  with  the  number  of  de¬ 
grees  of  freedom).  Exact  diagonalization  approaches 
(e.g.,  the  Lanczos  method)  suffer  from  this  exponen¬ 
tial  “catastrophe”.  Viewed  the  other  way  around,  this 
basis  set  scaling  is  what  restricts  today’s  classical  com¬ 
puter  to  simulating  only  small  quantum  computers. 
This  point  seems  indisputable,  but  should  not  be  taken 
as  proof  that  quantum  systems  cannot  be  simulated 
on  a  classical  computer.  By  the  same  token,  the  re¬ 
cent  claims  [4]  that  quantum  computers  can  simulate 
all  quantum  systems  efficiently  lacks  explicit  and  de¬ 
tailed  algorithms  for  specific  problems,  and  lacks  a 
generic  model  of  quantum  computation  including  the 
unitary  maps  (quantum  gates)  that  can  be  physically 
implementable.  Even  if  a  quantum  computer  existed, 
some  interesting  quantum  problems,  such  as  finding 
the  ground  state  of  a  general  quantum  Hamiltonian, 
do  not  yet  have  efficient  quantum  algorithms.  Finding 
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such  a  quantity  for  small  systems  is  relatively  routine 
on  a  classical  computer. 

Feynman  in  fact  analyzed  two  alternatives  for  sim¬ 
ulating  physics  with  computers  [3].  One  uses  a  proba¬ 
bilistic  classical  computer  that  would  produce  from  the 
same  input  as  given  to  a  physical  system  the  same  dis¬ 
tribution  of  outputs  as  observed  for  the  physical  sys¬ 
tem.  The  other  uses  a  computer  constructed  of  dis¬ 
tinctively  quantum  mechanical  elements  that  obey  the 
laws  of  quantum  mechanics.  This  latter  proposal  is  the 
quantum  computer. 

To  the  question,  “Can  quantum  systems  be  proba¬ 
bilistically  simulated  by  a  classical  computer?”,  Feyn¬ 
man’s  answer  was  unequivocally  “No”.  1 *  This  answer 
is  surprising  for  even  at  that  time  some  quantum  sys¬ 
tems  were  being  very  successfully  simulated  proba¬ 
bilistically  on  classical  computers,  mainly  by  quan¬ 
tum  Monte  Carlo  (QMC)  methods  [5].  To  the  question, 
“Can  quantum  systems  be  simulated  with  a  quantum 
computer?”,  his  answer  was  a  qualified  “Yes”.  He  be¬ 
lieved  almost  certainly  that  this  could  be  done  for  a 
system  of  bosons  but  was  unsure  that  it  could  be  done 
for  a  system  of  fermions.  In  this  paper  we  present  a  de¬ 
sign  for  a  universal  quantum  computer  that  will  sim¬ 
ulate  a  system  of  fermions.  Before  doing  so,  we  first 
discuss  some  problems  that  can  be  solved  by  a  proba¬ 
bilistic  simulation  of  a  quantum  system  on  a  classical 
computer  and  others  that  cannot. 

Probabilistic  simulations  of  quantum  systems  on 
a  classical  computer  are  mainly  performed  with  the 
use  of  the  Monte  Carlo  method.  These  statistical  ap¬ 
proaches  were  introduced  to  overcome  the  difficulty 
of  exponentially  growing  phase  spaces  by  numerically 
evaluating  the  accompanying  many-dimensional  inte¬ 
grals  by  sampling  from  a  function  assumed  to  be  non¬ 
negative.  On  a  classical  computer  one  can  probabilis¬ 
tically  simulate  a  quantum  system  like  liquid  He4  [6] 
and  produce  results  that  accurately  compare  with  ex¬ 
periment.  The  situation,  however,  is  far  from  satis¬ 
factory.  An  unsatisfactory  state  of  affairs  results  from 
the  frequent  breakdown  of  the  non-negativity  assump¬ 
tion  and  is  called  “the  sign  problem”.  The  sign  prob¬ 
lem  is  manifested  by  the  seemingly  exponentially  hard 
task  of  estimating  the  expectation  value  of  an  ob¬ 


1  There  is  as  yet  no  mathematical  proof  that  this  is  the  correct 

answer. 


servable  with  a  given  error.  Interestingly,  Feynman’s 
negativism  about  quantum  systems  being  probabilis¬ 
tically  simulated  by  classical  computers  was  a  claim 
that  negative  probabilities  were  unavoidable  because 
of  the  “hidden  variable”  problem  and  the  possible  vi¬ 
olation  of  Bell  inequalities.  The  extent  to  which  the 
sign  problem  is  a  hidden  variable  problem  is  unclear. 
On  the  other  hand,  QMC  methods  do  not  faithfully  ad¬ 
here  to  Feynman’s  idea  of  a  probabilistic  computer. 
Two  important  differences  are  that  most  QMC  simula¬ 
tions  are  non-local  and  performed  in  imaginary  time. 
Feynman  discussed  real-time  simulations  on  a  local 
computer.  Implications  of  these  differences  have  been 
noted  by  Ceperley  [7]  who  suggests  Feynman  really 
argues  just  against  simulating  quantum  dynamics  on  a 
local  classical  computer.  In  any  case,  the  known  prob¬ 
abilistic  simulations  on  a  classical  computer  clearly 
do  not  qualify  as  a  universally  efficient  computational 
scheme  for  general  quantum  many-body  problems. 
The  limiting  factors,  for  whatever  reasons,  arc  nega¬ 
tive  or  complex-valued  probabilities  whether  the  sim¬ 
ulations  are  done  in  real  or  imaginary  time. 

To  place  the  sign  problem  in  a  better  perspective, 
we  will  start  with  a  real-time  analysis  of  a  collection 
of  interacting  quantum  particles.  Quantum  mechanics 
tells  us  that  these  particles  either  obey  Bosonic  sta¬ 
tistics,  whereby  the  wave  function  is  symmetric  with 
respect  to  the  exchange  of  the  states  of  any  two  par¬ 
ticles,  or  obey  Fermionic  statistics,  whereby  the  wave 
function  is  antisymmetric  (changes  sign)  with  respect 
to  the  exchange  of  any  two  particles  [8].  Examples  of 
bosons  are  photons  and  gluons;  examples  of  fermi¬ 
ons  are  electrons,  protons,  neutrons,  and  quarks.  Of¬ 
ten  these  two  quantum  statistics  conveniently  and  ef¬ 
ficiently  map  onto  a  third,  quantum  spin  statistics  [9]. 
Still  in  other  cases,  when  particle  exchange  is  unlikely, 
particle  statistics  is  simply  ignored. 

For  a  given  initial  quantum  state  Itf'(O)),  a  quan¬ 
tum  computer  solves  the  time-dependent  Schrodinger 
equation 

d\V) 

\h^-  =  H\V)  (1) 

by  incrementally  propagating  the  initial  state  via 
|  </>  (r))  =  e_iA' w//l  . . .  Q->AiH/h  |  ^  (2) 

M  factors 

(t  =  M  At  and  the  Hamiltonian  H  is  assumed  time 
independent.)  It  should  be  reasonably  apparent  that  if 
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the  Monte  Carlo  method  is  applied  to  the  evaluation 
of  the  right-hand  side  of  this  equation,  it  is  faced 
with  sampling  from  oscillatory  integrands  that  are  not 
always  positive  and  have  unknown  nodal  surfaces. 
Further,  as  time  t  increases,  the  integrand  fluctuates 
with  increasing  rapidity.  While  clever  stationary-phase 
forms  of  the  QMC  method  have  been  developed, 
acceptable  solutions  are  possible  only  for  relatively 
short  times.  This  form  of  the  sign  problem  is  called  the 
dynamical  sign  problem ,  and  we  are  unaware  of  any 
efficient  [10]  real-time  QMC  simulations  for  bosonic, 
fermionic,  or  quantum  spin  systems. 

Years  ago  [5],  before  quantum  computers  were  pro¬ 
posed,  it  was  realized  that  by  transforming 
Schrodinger’s  equation  to  imaginary-time  r  via  f  — ► 
—\t\T  the  problem  with  the  rapid  fluctuations  was 
eliminated.  With  this  transformation,  called  Wick’s  ro¬ 
tation,  one  solves  the  diffusion-like  equation 


by  incrementally  propagating  the  initial  state  via 
|^(t))  =  e~Ar/y  •••e~AT//J^(0)).  (4) 

M  factors 

(r  =  M  Ar  and  the  Hamiltonian  H  is  assumed  time  in¬ 
dependent.)  This  transformation  permits  QMC  simula¬ 
tions  of  time-reversal  invariant  interacting  boson  sys¬ 
tems  to  a  high  degree  of  accuracy.  For  systems  of  in¬ 
teracting  quantum  spins  and  fermions  (or  bosons  with 
complex  Hermitian  Hamiltonians  [11]),  the  transfor¬ 
mation  does  not  solve  the  sign  problem.  For  quan¬ 
tum  spin  systems,  the  difficulty  is  finding  a  basis  in 
which  all  matrix  elements  of  the  positive-definite  op¬ 
erator  exp(— At H)  are  positive.  Most  notably  this 
difficulty  occurs  for  frustrated  quantum  spins.  For 
fermion  systems,  the  problem  is  the  Monte  Carlo 
process  causing  state  exchanges  that  because  of  the 
anti-symmetrization  requirement  just  happen  to  pro¬ 
duce  samples  which  are  as  frequently  positive  as  neg¬ 
ative.  For  the  sign  problem  found  in  both  types  of 
systems,  the  statistical  error  of  the  measured  observ¬ 
ables  grows  exponentially  fast  with  increasing  system 
size.  Another  difficulty  with  the  imaginary-time  ap¬ 
proach  is  analytically  continuing  the  results  back  to 
real-time  if  real-time,  i.e.  truly  dynamical,  information 
is  needed  [12].  This  continuation  is  an  ill-posed  prob¬ 
lem  whose  solution  places  extraordinary  demands  on 
the  simulation  [13]. 


In  this  paper,  we  will  focus  on  the  dynamical  sign 
problem  for  a  system  of  fermions,  seemingly  the  most 
challenging  case.  Eventually  we  will  give  a  detailed 
implementation  of  a  simulation  of  the  dynamical  prop¬ 
erties  of  a  collection  of  interacting  fermions  on  a  quan¬ 
tum  computer.  The  implementation  avoids  the  sign 
problem.  First,  in  Section  2  we  will  discuss  more  fully 
the  mathematical  origin  of  the  dynamical  sign  problem 
in  classical  computation  and  show  why  a  quantum  al¬ 
gorithm  overcomes  the  problem.  In  Section  3  we  will 
give  the  elements  required  for  Deutsch’s  quantum  net¬ 
work  model  of  a  quantum  computer  [14].  The  quantum 
gate  in  this  model  conveniently  allows  the  propagation 
of  systems  of  local  two  state  objects,  e.g.,  a  localized 
quantum  spin-^  particle  called  qubit.  We  also  propose 
a  universal  set  of  quantum  gates  (unitary  operators) 
that  allows  generic  propagation  of  systems  of  fermi¬ 
ons  (the  fabled  “Grassmann  chip”  [15]).  The  resulting 
fermion  algebra  has  been  the  main  technical  tool  for 
studying  the  classical  Ising  model  in  two  spatial  di¬ 
mensions  [16],  a  prototype  lattice  model  that  had  an 
enormous  impact  on  our  understanding  of  phase  tran¬ 
sitions.  Next,  in  Section  4,  we  show  how  this  propaga¬ 
tion  can  be  effected  by  the  quantum  spin  gate.  We  will 
demonstrate  the  polynomial  scaling  of  the  construc¬ 
tion  of  the  initial  state,  its  subsequent  time  propaga¬ 
tion,  and  the  measurement  of  some  observable.  Here 
we  will  also  demonstrate  the  control  of  the  error  in  the 
results.  In  Section  5,  we  apply  our  model  of  dynami¬ 
cal  fermion  computation  to  a  toy  problem  to  illustrate 
our  procedures  in  more  detail.  Finally,  in  Section  6, 
we  summarize  and  make  some  remarks  about  future 
research  directions. 

Our  universal  fermion  gate  and  its  mapping  to  the 
standard  universal  quantum  gate  is  similar  to  the  one 
recently  discussed  by  Bravyi  and  Kitaev  [17]  who 
actually  propose  that  a  quantum  computer  built  from 
fermions  might  be  more  efficient  than  one  built  from 
distinguishable  two  state  systems. 

2.  Dynamical  sign  problem 

In  order  to  understand  the  mathematical  origin  of 
the  dynamical  sign  problem  we  use  the  Feynman  path 
integral  formulation  [18]  for  continuum  systems  in 
the  first  quantization  representation.  In  this  formalism 
one  maps  a  quantum  problem  in  D  dimensions  into 
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a  classical  one  in  D  4-  1  dimensions  and  then  simu¬ 
lates  that  problem  probabilistically  on  a  classical  com¬ 
puter.  The  algorithm  is  efficient  except  for  the  repeti¬ 
tion  needed  to  obtain  sufficiently  good  statistics.  The 
“distinguishable  particle”  quantum  mechanical  prop¬ 
agator  of  a  system  represented  by  the  Hamiltonian 
H  =  ^  pf  +  V  (Tl)  is  expressed  as  [19] 

G(7l  ->  72/;  t)  =  (7?/,  f|e_i///|7£,  0) 

K(t)=n' 

=  J  V[R.(t)]cism,)],  (5) 

n(0)=n 

where  the  measure 

V[Jl(tj\  =  lim  Unit /M)~m 0/2 dTl\ 

L  J  M->oo 

and  the  action 

S[n(t)]  =  j drji(^)2- V(7l(t))j.  (6) 

0 

Bosonic  or  fermionic  statistics  are  introduced  by  ap¬ 
plying  the  corresponding  symmetrization  operator  to 
the  propagator,  Eq.  (5).  However,  because  the  dynam¬ 
ical  sign  problem  occurs  for  any  particle  statistics,  we 
will  ignore  particle  statistics  for  the  sake  of  simplic¬ 
ity. 

The  description  of  the  properties  of  different  phys¬ 
ical  systems  in  terms  of  correlations  of  physical  ob¬ 
servables  is  the  natural  way  to  compare  with  available 
experimental  information.  In  this  regard,  linear  re¬ 
sponse  theory  provides  a  way  to  compute  the  response 
of  a  system  to  a  weak  external  dynamical  perturba¬ 
tion  [20].  This  linear  response  is  always  expressed  in 
terms  of  a  time  correlation  function  of  the  dynamical 
variables  that  couple  to  the  perturbation.  For  exam¬ 
ple,  if  we  were  to  apply  an  external  time-dependent 
magnetic  field  and  we  wanted  to  calculate  the  average 
induced  magnetization,  we  would  have  to  compute  a 
time-dependent  magnetization-magnetization  correla¬ 
tion  function.  The  two-time  correlation  function  be¬ 
tween  arbitrary  local  dynamical  variables  A  and  B 
is 

CAB(t)  =  (A(t)B(O))  =  (eiw'  Ac~iH,B),  (7) 

if  the  Hamiltonian  is  time  independent.  Generically,  a 
stochastic  estimate  of  Cab(0  is 


Cab(0 


{A{-RM+x)B{TL\)e0^)P 


(8) 


where  the  configurations  {7 £/}  are  sampled  from  the 
probability  distribution  P  (positive  semidefinite  mea¬ 
sure),  and  <P  is  a  real-valued  function.  One  imme¬ 
diately  sees  that  the  origin  of  the  dynamical  sign 
problem  is  the  oscillatory  phase  factor  e10  that  leads 
to  large  phase  fluctuations  at  long  times.  Manifestly, 

| (Qi<P({7li }))p|  — ►  0  in  an  exponential  fashion  as  t  gets 
larger.  Therefore,  the  total  statistical  error  for  the  eval¬ 
uation  of  Cab(0  grows  exponentially  with  time  be¬ 
cause  of  large  cancellations  both  in  the  numerator  and 
denominator.  The  so-called  “fermion  sign  problem”  is 
a  particular  case  of  this  problem  when  e10  =  ±  1  and 
time  is  imaginary  [21]. 

Will  a  quantum  computer  solve  this  problem?  One 
often  hears  that  it  will  because  a  quantum  computer 
is  a  physical  system,  whether  a  system  of  fermions 
or  not,  and  physical  systems  have  no  dynamical  or 
fermion  sign  problems.  Furthermore  it  has  been  ar¬ 
gued  that  there  are  means  for  mapping  most  physi¬ 
cal  systems  to  a  quantum  computer  in  such  a  way  that 
the  quantum  computer’s  controlled  evolution  mimics 
that  of  the  physical  system  [3,22].  A  closer  look,  how¬ 
ever,  makes  the  situation  less  clear.  A  quantum  com¬ 
puter  is  a  computer,  and  as  such  it  suffers  from  lim¬ 
ited  accuracy.  More  importantly  this  type  of  computer 
predicts  results  stochastically,  meaning  each  measure¬ 
ment  is  one  member  of  the  ensemble  of  measurements 
possible  from  a  distribution  specified  by  the  modulus 
squared  of  the  wave  function  for  the  Hamiltonian  H 
modeled  by  the  quantum  computer.  For  a  fixed  physi¬ 
cal  time  t  >  0,  how  accurate  is  an  individual  measure¬ 
ment,  how  accurate  is  the  expectation  value  of  these 
measurements,  and  how  controlled  is  their  estimated 
variance?  Is  the  level  of  accuracy  and  control  achiev¬ 
able  polynomially  with  complexity  and  tl 

There  is  an  area  where  a  problem  similar  to  the  sign 
problem  has  been  recognized  and  resolved  by  quan¬ 
tum  computation.  Recently  it  was  shown  that  quan¬ 
tum  computation  is  polynomially  equivalent  to  clas¬ 
sical  probabilistic  computation  with  an  oracle  for  es¬ 
timating  the  value  of  simple  sums  of  rational  num¬ 
bers  called  quadratically  signed  weight  enumerators 
(QWGTs)  [23].  In  other  words,  if  these  sums  could  be 
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evaluated,  one  could  use  them  to  generate  the  quan¬ 
tum  statistics  needed  to  simulate  the  desired  quan¬ 
tum  system.  More  specifically,  what  was  demonstrated 
was  the  obtainability  of  expectation  value  of  opera¬ 
tors  in  quantum  computation  by  evaluating  sums  of  the 
form 

S(A,B.x,y)=  (-l)bJBbxwyH~w,  (9) 

h:Ah= 0 

where  A  and  B  are  0-1 -matrices  with  B  of  dimen¬ 
sion  n  x  n  and  A  of  dimension  m  x  n.  The  vari¬ 
able  b  in  the  summand  ranges  over  0-1 -column  vec¬ 
tors  of  dimension  n,  bT  denotes  the  transpose  of 
b ,  \b\  is  the  weight  of  b  (the  number  of  ones  in 
the  vector  b),  and  all  calculations  involving  A ,  B 
and  b  are  done  modulo  2.  The  absolute  value  of 
S(A,  B,xyy)  is  bounded  by  (\x\  +  \y\)n.  Quantum 
computation  corresponds  to  the  problem  of  determin¬ 
ing  the  sign  of  S (A ,  \t(A) ,  k  J)  with  the  restrictions 
of  having  dg(A)  =  /,  k  and  /  being  positive  inte¬ 
gers,  and  |S(A,  lt(A),&,/)|  ^  ( k 2  +  l2)n/2/2.  dg(A) 
is  a  diagonal  matrix  formed  from  the  diagonal  el¬ 
ements  of  A  and  lt(A)  is  a  lower  triagonal  ma¬ 
trix  formed  from  the  lower  triangular  elements  of 
A.  Details  of  this  quantum  algorithm  can  be  found 
in  [23]. 

The  main  point  is  that  these  sums  are  similar  to 
the  numerator  of  Eq.  (8),  and  attempts  to  estimate 
them  by  random  sampling  result  in  exponentially  bad 
signal  to  noise  ratios.  In  the  case  of  QWGTs,  quantum 
computers  can  estimate  the  sum  exponentially  better 
than  classical  computers,  but  the  estimate  is  not 
exact.  The  situation  for  the  dynamical  sign  problem 
is  similar:  Quantum  computers  cannot  obtain  exact 
values  for  the  desired  correlation  functions,  but  can 
obtain  estimates  sufficiently  exact  to  avoid  the  sign 
problem  suffered  by  the  known  classical  algorithms 
and  to  yield  usable  information  about  the  physical 
models  simulated. 

In  this  paper  we  will  show  explicitly  how  the  sign 
problem  is  avoided  in  the  case  of  simulating  fermi¬ 
ons.  Below  we  will  give  a  means  for  translating  lo¬ 
cal  fermion  Hamiltonians  into  the  Hamiltonians  avail¬ 
able  in  the  standard  model  of  quantum  computation. 
In  contrast  to  quantum  simulations  on  a  classical  com¬ 
puter  this  translation  prevents  uncontrolled  exchange 
processes  that  are  the  dominant  source  of  the  fermion 
sign  problem.  With  respect  to  the  dynamical  sign  prob¬ 


lem,  we  then  argue  by  using  standard  error  correction 
analysis  developed  for  the  standard  model  of  quantum 
computing  that  these  gates  will  enable  sufficiently  ac¬ 
curate  measurements  of  correlation  functions  so  the 
accuracy  of  the  average  of  these  measurements  will 
be  dominated  by  the  statistical  error.  The  statistical  er¬ 
ror  is  problem  dependent  but  polynomially  bounded, 
so  that  the  difficulty  associated  with  phase-weighted 
averages  is  eliminated. 

3.  Models  of  quantum  computation 

The  quantum  control  model  of  quantum  computa¬ 
tion  assumes  the  existence  of  physical  systems  that 
can  be  controlled  by  modulating  the  parameters  of  the 
system’s  Hamiltonian  Hp.  The  control  possibilities 
are  abstracted  and  used  to  implement  specific  quan¬ 
tum  gates  that  represent  the  unitary  evolution  of  the 
physical  system  over  a  time  step  obtained  by  specific 
modulations  of  the  Hamiltonian.  In  most  treatments, 
the  physical  systems,  together  with  the  gates,  are  then 
taken  as  the  abstract  model  of  quantum  computation. 
The  quantum  control  and  quantum  gate  viewpoints 
are  effectively  equivalent,  but  to  tie  the  computational 
model  to  the  physics  simulation  problem  more  closely, 
we  choose  to  describe  quantum  computation  from  the 
point  of  view  of  quantum  control;  that  is,  we  will 
assume  an  Hp.  In  this  context  we  begin  by  giving 
the  standard  model  of  quantum  computation  and  then 
defining  an  alternative  model  based  on  fermions. 

Defining  a  model  of  quantum  computation  consists 
of  giving  an  algebra  of  operators,  a  set  of  controllable 
Hamiltonians  (Hermitian  operators  in  the  algebra),  a 
set  of  measurable  observables,  and  an  initial  state  of 
the  physical  system.  In  the  simplest  case,  the  observ¬ 
ables  are  measured  by  the  method  of  projective  mea¬ 
surements,  and  the  initial  state  of  the  physical  system 
is  an  expectation  value  of  the  algebra’s  operators. 

3. 1.  Standard  model  of  quantum  computation 

The  standard  model  of  quantum  computation 
(Deutsch’s  quantum  network  representation)  is  based 
on  an  assembly  of  two  state  systems  called  qubits ,  con¬ 
trolled  by  one-  and  two-qubit  Hamiltonians,  and  on 
a  measurement  process  determined  by  one-qubit  ob¬ 
servables. 
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Operator  algebra.  It  is  convenient  to  define  the 
standard  model  through  the  algebra  of  operators  acting 
on  the  qubits.  This  algebra  is  generated  by  the  unit  and 
Pauli  matrices  ox,  oy  and  crz  for  each  qubit  j. 


These  matrices  represent  quantum  operators  with 
mixed  commutation  relations  and  span  the  space  of 
complex-valued  2x2  matrices.  For  qubits  j  /  k ,  the 
o' s  commute,  and  for  qubits  j  =  k ,  they  satisfy  the 
relation  olJLov  -f  crvotl  =  2 8flvt  (/x,  v  =  xy  y,  z).  For  a 
quantum  register  with  n  qubits,  one  may  take  the  op¬ 
erator  ojL  defined  in  terms  of  a  Kronecker  product 

<7^  =  1  ®  1  ®  •  •  •  ®  o[L  ®  •  •  •  ®  l 
jth  factor 

of  matrices  acting  on  n  two-dimensional  linear  spaces. 
Thus  ofi  admits  a  matrix  representation  of  dimension 
2n  x  2n. 

Control  Hamiltonians.  Control  of  qubits  is  achieved 
by  applying  Hamiltonians  that  act  on  either  one  or 
two  qubits.  A  theorem  [24,25]  in  quantum  information 
processing  says  that  a  generic  operation  on  a  single 
qubit  and  any  interaction  between  two  qubits  is  suffi¬ 
cient  for  building  any  unitary  operation.  We  take 

Hp(t)  =  £>,  +  aVj  (0^  ] 

j 

+  J2a‘j  ^azaZ' 

ij 

where  the  afl(t)  are  controllable.  Ideally,  no  con¬ 
straints  on  the  control  functions  are  assumed.  How¬ 
ever,  it  is  often  simpler  to  design  the  required  con¬ 
trol  by  assuming  that  only  one  of  the  a, ft)  is  non¬ 
zero  at  any  time.  A  quantum  algorithm  for  this 
model  of  quantum  computation  consists  of  prescrib¬ 
ing  the  control  functions  [26].  A  convenient  measure 
of  the  complexity  of  such  an  algorithm  is  the  integral 

Jo  d (^e  action  of  the  algorithm).  The 
quantum  gates  are  simply  specific  unitary  evolutions 
that  may  be  implemented  in  terms  of  Hp.  A  conve¬ 
nient  universal  set  of  gates  is  given  by  operators  of  the 


form  exp (itr^Tr/4)  and  exp(io\!crJ;r/8).  In  the  quan¬ 
tum  network  representation  of  the  standard  model,  an 
algorithm  is  a  specific  sequence  of  these  operators  ap¬ 
plied  to  the  initial  state  of  the  qubits. 

Initial  state.  The  initial  state  of  the  qubits  is  assumed 
to  be  an  n  term  Kronecker  product  of  the  state  |0)  = 
(q)  which  is  an  eigenstate  of  oz  with  eigenvalue  1. 
The  state  is  completely  determined  by  the  expectation 
values  (Ola^.lO),  which  arc  1  if  the  olflj  are  all  o[  or 
the  identity,  and  are  0  otherwise.  Physically,  the  initial 
state  has  all  “spins”  up. 

Measurement.  The  final  feature  of  the  model  of  com¬ 
putation  is  the  specific  means  for  extracting  informa¬ 
tion  after  a  sequence  of  operations  has  been  applied  to 
the  initial  state.  In  the  standard  model,  it  is  always  pos¬ 
sible  to  apply  a  projective  (von  Neumann)  measure¬ 
ment  [27]  using  the  observables  of  With  this  capa¬ 
bility,  it  is  unnecessary  to  give  an  initial  state  explic¬ 
itly,  as  the  desired  state  can  be  prepared  by  using  mea¬ 
surement  and  operations.  To  learn  the  expectation  of 
an  observable  at  the  end  of  an  algorithm,  one  repeats 
the  algorithm  and  measurement  procedure  many  times 
and  averages  over  the  measurements  until  the  desired 
accuracy  is  achieved. 

For  a  description  of  the  standard  model  of  quantum 
computation  in  terms  of  quantum  Turing  machines, 
see  [28].  Quantum  networks  are  discussed  in  [24]. 
Introductory  descriptions  of  the  standard  model  may 
be  found  in  [29,30]. 

3.2.  Fermion  model  of  quantum  computation: 
Grassmann  chip 

Somewhat  analogously,  we  now  describe  a  standard 
model  of  fermion  computation.  For  simplicity  wc 
only  consider  spinless  fermions,  i.e.  fermions  without 
internal  spin  degrees  of  freedom,  although  we  could 
have  considered  more  general  fermionic  algebras  with 
internal  degrees  of  freedom  [9].  Physically,  a  system  of 
spinless  fermion  might  be  a  system  of  spin-  \  electrons 
in  a  magnetic  field  sufficiently  strong  to  polarize  it 
fully.  The  basic  system  of  this  model  is  a  state  (or 
fermionic  mode)  that  can  be  occupied  by  0  or  1 
spinless  fermion.  We  define  the  model  for  n  such 
modes. 
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Operator  algebra.  We  define  the  model  through  the 
algebra  of  the  spinless  fermion  operators  aj  and  aj  for 
each  qubit  j  (7  =  1, i.e.  through  the  algebra 
of  2 n  elements  satisfying  canonical  anticommutation 
relations 

[aitaj)  =  0,  =8ijy 

where  [A,  B)  =  AB  +  BA  denotes  the  anticommuta¬ 
tor  or  Jordan  product,  a J  (aj)  creates  (annihilates)  a 
spinless  fermion  in  state  (mode)  j.  Each  element  ad¬ 
mits  a  matrix  representation  of  dimension  2n  x  2n .  The 
fermion  algebra  is  isomorphic  (as  a  *-algebra)  to  the 
standard  model  (or  Pauli)  algebra.  The  isomorphism  is 
established  through  the  Jordan-Wigner  mapping  [31]. 

Control  Hamiltonians.  We  take 

HP  =  J2[<*j(t)aj  +6tj(t)a]] 

j 

+  52[a<7(0(a,V/  +a]ai)  +  Pij(t)a]aia]aj]. 
ij 

This  is  a  universal  Hamiltonian,  i.e.  any  other  Hamil¬ 
tonian  for  a  system  of  interacting  spinless  fermions 
can  be  generated  by  it.  Physical  operators  must  be 
(Hermitian)  products  of  even  degree  involving  combi¬ 
nations  of  the  creation  and  annihilation  operators  such 
as  the  terms  in  the  last  two  summands  of  the  Hamil¬ 
tonian  above. 

Initial  state.  The  initial  state  is  assumed  to  be  an  n 
term  Kronecker  product  of  the  state  |0)  which  is  an 
eigenstate  of  the  number  operator  a^aj  with  eigen¬ 
value  0.  The  state  is  completely  determined  by  the  ex¬ 
pectation  values  (0|tftfl7|0)  =  0  for  all  j.  Physically, 
the  initial  state  has  all  modes  unoccupied. 

Measurement.  Measurements  can  again  be  performed 
by  using  von  Neumann’s  scheme  of  projective  mea¬ 
surements.  In  Section  4.3,  we  will  discuss  another 
scheme  more  appropriate  for  the  physical  systems  at 
hand. 

In  the  next  subsection  we  show  how  to  simulate  the 
fermion  model  by  using  the  standard  spin-j  model. 
In  particular  it  is  possible  to  efficiently  map  the 
fermion  Hamiltonians  to  Pauli  operators  which  can 
be  simulated  using  the  control  Hamiltonians  of  the 


standard  model.  This  establishes  that  these  two  models 
of  computation  are  polynomially  equivalent.  Here  the 
point  of  view  is  similar  to  the  one  used  for  classical 
models  of  computation:  the  simulation  of  one  model 
by  another  establishes  their  equivalence. 


4.  Fermion  computation  via  the  standard  model 

In  the  previous  Section  we  gave  the  elements  re¬ 
quired  for  Deutsch’s  quantum  network  model  of  a 
quantum  computer  [  1 4]  and  proposed  a  universal  set  of 
quantum  gates  (unitary  operators)  that  allows  generic 
propagation  of  systems  of  fermions  (the  fabled  “Grass- 
mann  chip”  [15]).  Here  we  show  how  this  propaga¬ 
tion  can  be  effected  by  the  quantum  spin  gate.  We  will 
demonstrate  the  polynomial  scaling  of  the  construc¬ 
tion  of  the  initial  state,  its  subsequent  time  propaga¬ 
tion,  and  the  measurement  of  some  observable.  We 
will  also  demonstrate  the  control  of  the  error  in  the 
results. 

The  first  step  is  the  observation  that  the  set  of  2 n 
matrices  y ^  (of  dimension  2n  x  2n )  satisfying  the 
Clifford  algebra  identities 

{Yp>  Yv)  =  28tiv 

admits  a  representation  in  terms  of  Pauli  matrices 
(Brauer-Weyl  construction) 

Y\  =  <*x  -  Y2  =  . 

Yi  =°'ax<  Y4  =  cr^Oy, 

Y5  =  a}aj ,  Y6  =  crJ  a}a^ , 


Yin- 1  = 


n  —  1 

FK 

L  y=l 


/j 

°x  ’  Y2n  — 


n- 1 

FK 

L  ;=1 


V 


The  following  mapping  of  fermion  operators 


fj- 1 


ca 


f[ -CTz  F-  =  (“O'  x  a\o\  —  al  1 a‘_ 


,1  =  1 


=  (-iy 


-1  Y2j-\  -lY2j 


(-iy  yv.2 ■  ■  t/  '(t2 
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where  oJ±  =  (gx  ±  \Oy)/ 2,  defines  a  *-algebra  iso¬ 
morphism  to  the  algebra  of  operators  of  the  standard 
model.  It  is  the  so-called  spin-|  Jordan-Wigner  trans¬ 
formation  [31],  and  has  the  property  that  hj  =  a^.aj  = 

g\_gj_  =  j(l  +  g[).  We  note  that  hj  is  a  “local”  par¬ 
ticle  number  (or  density)  operator  and  many  types  of 
interaction  in  physical  systems  are  of  the  form  “den¬ 
sity  times  density”  which  simplifies  the  simulation  as 
we  will  see. 

It  is  important  to  emphasize  that  the  success  of 
our  approach  depends  upon  the  mapping  of  algebras 
(and  not  of  Hilbert  spaces).  In  this  regard  it  is  relevant 
to  mention  that  the  transformation  just  presented  is 
a  particular  case  of  a  more  general  set  of  mappings 
that  we  would  like  to  name  generalized  Jordan- 
Wigner  transformations  [9].  It  is  possible  to  imagine 
a  quantum  computer  implemented,  for  example,  with 
a  three  state  unit  (5=1)  instead  of  a  qubit.  In  such 
a  case,  these  generalized  transformations  still  allow 
one  to  simulate  fermions  or  particles  with  arbitrary 
statistics. 

Two  additional  comments  are  in  order:  The  map¬ 
ping  for  aj  and  aj  described  above  corresponds  to 
a  one-dimensional  array  of  spins.  The  extension  to 
higher  spatial  dimensions  can  be  done  [9,32,33]  in  var¬ 
ious  ways.  A  straightforward  extension  to  two  dimen¬ 
sions  is  to  re-map  the  sites  of  a  two-dimensional  array 
onto  a  one-dimensional  string  and  proceed  as  before. 
Also  there  is  nothing  special  about  using  the  fermion 
instead  of  a  quantum  spin  as  an  alternative  model  of 
computation.  One  could  have  just  as  well  used  the  hard 
core  boson  [9].  The  main  question  is  whether  different 
algebras  admit  a  physical  realization.  For  hard-core 
bosons  this  realization  is  He4  atoms. 

4.1.  Evolution 

Given  a  fermion  model  algorithm,  it  is  necessary 
to  efficiently  obtain  a  corresponding  standard  model 
algorithm  that  at  least  approximates  the  desired  evo¬ 
lution.  The  general  principle  is  to  map  the  time- 
dependent  fermion  Hamiltonian  H(t)  =  //,  to  the 

standard  model  operators  via  the  Jordan-Wigner  trans¬ 
formation,  express  the  result  in  terms  of  a  sum  of  sim¬ 
ple  products  of  Pauli  operators,  and  then  use  the  Trot¬ 
ter  approximation 


e-iA t(H0+Hi+  -)/h  _  +C>((A/)2).  (11) 

i 

Each  time  step  At  is  chosen  so  that  the  final  error 
of  the  simulation  is  sufficiently  small.  Provided  that 
the  number  of  terms  in  the  sum  is  polynomially 
bounded  in  the  number  n  of  qubits  or  fermionic  modes 
and  provided  that  each  term  can  be  polynomially 
simulated,  the  simulation  is  efficient  in  n  and  1  /error. 

To  see  how  to  do  the  simulation,  consider  the 
example  of  the  bilinear  operator  Hc  =  a\a^  +  djci\  in 
the  control  Hamiltonian  of  the  fermion  model: 


He  =  ( - 1 );  [ct!  CT,1  ■■■ol  1 GJ 


,  _  1  _  1  _2  /'-l  /I 

-}-  G z  G+Gz  •  •  •  Gz  G _  J 


(-n7'i 


r  12  /'-i  y  ,  12  ;-i  n 

[Gx  Gz  •  •  •  Gz  Gx  4-  GyGz  •  •  •  G:  Gy  J 


J-lJl 


It  is  readily  checked  that  the  Jordan-Wigner  transfor¬ 
mation  for  the  other  terms  in  the  control  Hamiltonians 
are  also  decomposable  into  sums  of  a  few  products  of 
Pauli  operators. 

The  whole  idea  of  a  quantum  computer  is  simu¬ 
lating  the  operations  we  want  by  using  unitary  ma¬ 
trices  U  =  exp(— iA tHp/h).  These  unitary  matrices, 
representing  quantum  gates,  perform  reversible  com¬ 
putation  and  are  case-dependent.  For  our  particular 
case,  we  know  how  to  simulate  H  =  gz  in  the  spin-^ 
case  (it  is  directly  implemented  in  the  standard  model), 
so  we  ask  what  set  of  unitary  operations  produce 
the  evolution  U  =  exp(— iA tHJh).  In  other  words, 
how  do  we  write  a  U  =  U\  . . .  £/*  such  that  Hc  = 

HU?  Consider  for  example  the  Hamiltonian  Hx  = 
g\gz  •  •  •  o’/1  gx  .  The  procedure  is  as  follows:  The 
unitary  operator 


Ux  =c'W  =  j^[i  +  ^] 


1 

Ti 


(-,  0 


®  1  ®  ®  1 


takes  a)  — ►  a\ ,  i.e.  U\a}U \  —  <Tt' .  The  operator 

U2  =  e1?^2  =  2=1!  +  icr’cr2] 

V2 

takes  g\  — ►  cr^cr?.  The  next  step  is 


C/3=eiW 
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to  take  Qy oj  -*  —g\g^g}.  By  successively  similar 
steps  we  easily  build  the  required  string  of  operators: 
a l a}  ■  ■  •  <r/- ' cr/ . 

If  j  is  odd, 

Uj=eiiaM 

will  take  a' <72  •  •  •  al ~ 1  (-l)l^_1)/2|o-Jcrz2  •  •  -cr/, 

where  [m/l]  is  the  integer  part  of  m/l.  The  final 
operator 

Uj+ 1  =eJW 

will  bring  the  control  operator  to  the  desired  one  (up 
to  a  global  phase  (— l)f0“l)/2]): 


If  j  is  even,  we  need  an  additional  unitary  operator  that 
flips  the  first  qubit’s  crvl  into  a  aj .  This  flip  is  achieved 
with  the  operator 

Uj+  2  =  o  ^  ®  1  ®  ®  1 

that  takes  Gy  ->  crj . 

Hence,  to  construct  this  non-local  fermion  operator 
from  the  standard  model  requires  additional  steps  that 
are  proportional  to  j .  This  number  scales  polynomi¬ 
al^  with  the  complexity  so  the  construction  is  efficient 
if  the  standard  model  is  efficient. 

The  one  and  two-body  nature  of  naturally  occurring 
interactions  means  that  a  term  in  a  second-quantized 
representation  of  a  Hamiltonian  only  has  one  of  two 
forms:  either  ajaj  or  ajajala^  We  just  demonstrated 
how  to  handle  the  first  case.  The  second  case  merely 
requires  applying  that  algorithm  twice.  This  squares 
the  complexity. 

4.2.  State  preparation 

In  this  section  we  discuss  the  preparation  of  states 
of  physical  relevance.  Clearly,  the  preparation  of  the 
initial  state  is  a  very  important  step  since  the  study  and 
efficiency  of  the  given  physical  process  one  wants  to 
simulate  depends  upon  it. 

Consider  a  system  of  Ne  fermions  and  n  operators 

(single  particle  states).  A  generic  A^-particle  state 
of  a  Hilbert  space  'H^e  of  antisymmetrized  wave 


functions  can  always  be  expanded  in  terms  of  the 
antisymmetric  states 

Ne 

7  =  1 

where  fct  creates  a  state  j  and  |vac)  =  |0)  0  |0)  •  •  •  ® 
|0)  is  the  vacuum  state  (i.e.  bj\0)  =  0,  Vy).  The 
operator  is  in  general  a  linear  combination  of  a^’s, 
i.e.  b^  =  YHi=\a]Pij  where  P/;  is  some  matrix  and 
Ne  ^  n. 

The  states  \<Pa)  (a  =  1, . . . ,  ))  in  general  form 

an  overcomplete  set  of  non-orthogonal  states  that  span 
the  whole  HNe ,  i.e.  redundantly  generate  HNe  •  They 
are  known  as  Slater  determinants  [20].  Typically,  \<Pa) 
is  the  result  of  a  self-consistent  mean-field  (or  gen¬ 
eralized  Hartree-Fock)  calculation.  Even  a  Bardeen- 
Cooper-Schrieffer  superconducting  state,  which  does 
not  preserve  the  number  of  particles,  can  be  written  in 
this  way  after  an  appropriate  canonical  transformation 
which  redefines  the  vacuum  state  [34]. 

One  can  easily  prepare  the  states  \<Pa)  noticing  that 
the  quantum  gate,  represented  by  the  unitary  operator 

Um  = 

when  acting  on  the  vacuum  state,  produces  b]n  |0)  up 
to  a  phase  e1? .  Therefore,  the  successive  application  of 
similar  unitary  operators  will  generate  the  state  \<Pa) 
up  to  a  global  phase. 

Except  for  very  small  systems  the  total  Hilbert 
space  is  too  large  to  be  fully  used  (it  has  an  exponential 
growth  with  increasing  system  size).  In  practice,  one 
works  in  a  subspace  of  that  closely  represents  the 
physical  state  one  is  trying  to  simulate.  Generically,  as 
initial  state,  we  will  consider  a  very  general  expression 
of  a  many-fermion  state: 

N 

|  &  (t  =  0))  =  ^  ^  3 a  |  ^ or  )> 

or = 1 

where  the  integer  /V  is  a  finite  and  small  number.  The 
state  can  be  prepared  efficiently  (in  N)  by  a  number  of 
procedures.  We  now  describe  one. 

To  make  the  description  simple,  we  will  assume 
J2a=\  la<*l2  =  1  and  (&a\&fi)  =  <w,  which  is  equiv¬ 
alent  to  requiring  {|<£a)}  to  be  an  orthonormal  set  and 
|^(r  =  0))  to  be  normalized  to  unity.  With  these  as¬ 
sumptions  the  steps  of  the  state  preparation  algorithm 
are: 
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(1)  Adjoin  A  auxiliary  (ancilla)  qubits,  each  in  the 
state  |0),  to  the  vacuum  of  the  physical  system. 
The  resulting  state  is 

|0)  0  |0)  0  •  •  •  |0>  0|vac)  =  1 0} a  0  | vac).  (12) 

i  -  -  - 

N 

(2)  From  this  state  generate  ^a=\a<x\a)  ®  lvac) 
where  |a)  is  an  ancilla  state  with  only  the  a 
qubit  being  1 1 ).  The  procedure  for  generating  this 
combination  of  states  is  described  below. 

(3)  For  each  a  =  1 , . . . ,  A,  conditional  on  the  a  qubit 
being  |1),  apply  the  state  preparation  procedure 
for  | <Pa).  The  resulting  state  is 

N 

£a«|a>®|*«>.  (13) 

Ot=  1 

(4)  From  this  state  generate 
a<y|0)a  0  \<Pa)  4-  terms  without  |0)a. 

(14) 

This  step  will  also  be  described  below. 

The  final  step  is  accepted  if  a  measurement  shows  all 
the  ancillas  being  returned  to  |0)a.  The  probability  of 
successful  preparation  is  thus  5Za=i  |a»l2/A  =  1/A. 
Consequently,  on  average,  A  trials  will  be  needed 
before  a  successful  state  preparation. 

The  procedure  to  produce  step  2  is  most  easily 
described  by  example.  We  will  assume  A  =  2.  The 
problem  then  is  to  generate  a\  |  10)  ®  |vac)  +  a2 10 1 )  0 
| vac)  from  |00)  0  |vac).  In  what  follows  all  operations 
will  be  only  on  the  ancilla  part  of  the  initial  state  so  we 
will  not  explicitly  show  the  vacuum.  We  also  note  that 
one  can  always  apply  a  rotation  to  a  given  qubit  that 
will  take  |0)  into  jc |0)  4  y\  1)  with  |jc|2  +  |y|2  =  1.  The 
steps  of  the  procedure  are: 

(2.1)  Adjoin  an  ancilla  qubit  |b)  initially  being  |0). 
The  initial  state  is  now  |0)  0  |00). 

(2.2)  Conditional  on  |b)  =  |0),  rotate  the  or  =  1  qubit, 
and  then  conditional  on  the  ct  =  1  qubit  being 
|1),  flip  |b>: 

*i  |0)  0  |00)  +  yi  1 1)  0  1 10).  (15) 

(2.3)  Conditional  on  |b)  being  |0),  rotate  the  a  =  2 
qubit,  and  then  conditional  on  the  a  =  2  qubit 
being  |1),  flip  |b): 


^2|O)0|OO)+jciy2|l)®|Ol)  +  yi|l)0|lO). 

(16) 

(2.4)  Project  out  the  states  with  |b)  being  |1): 

*iwl01)  +  yi|10).  (17) 

The  rotations  are  chosen  so  that  a\  =  y \  and  a2  = 

x\y2- 

For  the  explanation  of  step  4,  we  will  display  the 
physical  states.  The  problem  is:  Starting  with  a\  |  10)  0 
|<Z>i)  4-  a2 101 )  0  |<£2),  produce  (14). 

(4.1)  Adjoin  an  ancilla  qubit  |b)  initially  being  |0). 
The  initial  state  is  now  a 1 10)  0  1 1 0)  0  |<Z>i)  4- 
a2|O)0  |O1)0|#2). 

(4.2)  Conditional  on  |b)  being  |0),  rotate  the  a  =  1 
qubit,  and  then  conditional  on  the  a  =  1  qubit 
being  |1),  flip  |b): 

a1(jt||0)®|00)  +  y1|l)®|10))®|<P|) 

+  a2(jri|0)®|01)  +  >i|l)®|ll>)®|tf»2>. 

(18) 

(4.3)  Conditional  on  |b)  being  |0),  rotate  the  a  =  2 
qubit,  and  then  conditional  on  the  a  =  2  qubit 
being  1 1),  flip  |b): 

aiJti(jt2|0)®  |00)  +  y2|l)®  |01>)  ®  |<f>j) 

+  32A-i(a:2|0)  ®  100)  +  J2l  1 )  ®  |01»  ®  |<f>2>- 

(19) 

(4.4)  Project  out  the  states  with  |b)  =  |0): 

x\X2(a\ |00)  0  |0i)  +  a2|00)  0  \4>2))-  (20) 

The  rotations  are  chosen  so  that  jci*2  equals  \/y/~N 
where  N  =  2.  Comparing  step  2  with  step  4,  one 
sees  they  are  structurally  identical,  differing  by  the  set 
of  amplitudes  generated  and  the  complementarity  of 
the  subspaces  selected  for  the  final  result.  This  latter 
difference  in  some  sense  makes  one  procedure  the 
inverse  of  the  other.  For  the  case  of  N  >  2,  one  simply 
replaces  steps  2.2  and  2.3  and  steps  4.2  and  4.3  by  “do 
loops”  over  a  from  1  to  A. 

On  average,  the  entire  procedure  needs  A  trials 
before  a  successful  state  preparation.  (In  many  cases, 
the  other  measurement  outcomes  can  be  used  also 
to  avoid  too  many  trials.)  Construction  of  the  initial 
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state  thus  scales  as  0(N2nNe)  ^  0(N2n2)  so  unless 
the  number  of  Slater  determinants  is  exponentially 
large,  general  many-fermion  states  can  be  initialized 
efficiently. 

4.3.  Measurement 

While  there  is  a  variety  of  physical  observables  one 
measures  experimentally  and  calculates  theoretically, 
at  this  time  it  is  difficult  to  demonstrate  that  they  all 
can  be  computed  efficiently  on  a  quantum  computer. 
Fortunately,  we  will  now  argue  that  one  important 
class  of  observables,  the  temporal  correlation  func¬ 
tions  Cab{0 ,  can  be  computed  not  only  efficiently  but 
also  accurately.  These  functions  describe  the  tempo¬ 
ral  evolution  of  some  observable  A(t)  in  response  to 
some  weak  external  stimulus  that  couples  to  the  sys¬ 
tem’s  variable  B( 0).  They  are  at  the  heart  of  under¬ 
standing,  for  example,  the  optical  properties  of  mate¬ 
rials. 

The  goal  is  to  determine  correlations  of  the  form 
CAB(t)  =  (A(t)B(O))  =  (e'Ht Ae~]Ht B)  up  to  a  suf¬ 
ficiently  small  statistical  error.  Clearly,  measuring  ef¬ 
ficiently  Cab  is  not  possible  for  an  arbitrary  A  and 
B.  One  sufficient  condition  is  that  A  and  B  are  effi¬ 
ciently  simulatable  Hamiltonians.  This  observation  is 
based  on  a  method  for  determining  Cab  refined  by  Ki- 
taev  [35]  and  applied  to  the  measurement  of  correla¬ 
tion  functions  by  Terhal  and  DiVincenzo  [36].  Here  we 
give  a  different  method  based  on  an  idea  given  in  [37]. 

A  general  principle  that  can  be  used  to  obtain 
Cab  is  to  decompose  the  operator  whose  expectation 
needs  to  be  determined,  i.e.  A(t)B{ 0),  into  a  small 
sum  of  operators  of  a  simpler  form  and  measure  each 
summand  individually.  Our  method  directly  measures 
expectations  of  the  form  {U^V)  when  algorithms 
for  implementing  the  unitary  operators  U  and  V 
are  available.  General  correlation  functions  are  then 
determined  by  decomposing  operators  using  a  unitary 
operator  basis,  for  example  the  one  consisting  of 
products  of  Pauli  matrices. 

The  method  for  measuring  {U^V)  consists  of  the 
following  steps: 

(1)  Adjoin  via  a  direct  product  an  ancilla  (i.e.  an 
auxiliary)  qubit  a  in  the  state  (|0)  +  1 1 ))/ \Jl  with 
density  matrix  pa  =  (1  +  erf)/ 2  to  the  initial  state 
of  the  system  described  by  the  density  matrix  p. 


(2)  Apply  the  conditional  evolutions  U\  =  |0)(0|  8> 
£/  +  |1>(1|  ®  1  and  U2  =  |1)(1|  ®  V  +  |0)(0|  8)  11 
( U  =  U\U2).  The  methods  of  [24]  may  be  used  to 
implement  these  evolutions  given  algorithms  for 
U  and  V. 

(3)  Measure  2<rf  =  erf  +  ierf  =  2|0)(1|.  This  may  be 
done  by  measuring  erf  and  erf  in  sufficiently  many 
independent  trials  of  these  steps. 

(4)  Given  the  initial  density  matrix  p,  the  expectation 

[°x  +  i<Ty)p,®p  =  2Tr„+i  [f/^|0)(l  |t7pa  ®  p] 

=  Tr„+1[|0><l|®f/tVAl®p] 

=  Tr n[CVp)  =  (Cv)p,  (21) 

as  desired.  The  statistical  noise  in  the  measure- 
ment  of  (U^V)P  is  determined  by  that  of  two  bi¬ 
nary  random  variables  and  therefore  depends  only 
on  the  value  of  Tr n[U^Vp],  which  is  inside  the 
unit  complex  circle.  As  a  result  it  is  a  simple  mat¬ 
ter  to  determine  the  number  of  measurement  at¬ 
tempts  required  to  achieve  sufficient  statistical  ac¬ 
curacy. 

The  procedure  for  measuring  Cab(0  can  now  be 
summarized  as  follows:  First  express  A  =  A(0)  and 
B  =  B{  0)  as  a  sum  of  unitary  operators  A  =  \  Aj 

and  B  =  Yl"yB=\  #/'•  A  convenient  unitary  operator  ba¬ 
sis  that  worxs  well  for  the  local  observables  of  interest 
consists  of  all  the  products  of  Pauli  operators,  as  each 
such  product  is  easily  implemented  as  a  quantum  al¬ 
gorithm.  Then,  for  each  j  and  j'  one  uses  the  just  de¬ 
scribed  method  with  U  =  e]Ht  A^Q~]Ht  and  V  =  Bj>  to 
obtain  {Aj(t)By( 0)).  V  may  be  implemented  by  sim¬ 
ulating  the  evolution  under  //,  applying  By,  and  then 
undoing  the  evolution  under  H . 

An  alternative  approach  to  the  measurement  process 
is  von  Neumann’s  projection  method.  We  sketch  it 
here  for  completeness  and  comparison.  In  this  ap¬ 
proach  we  also  add  an  auxiliary  (ancilla)  degree  of 
freedom  to  the  problem.  Suppose  that  this  extra  qubit 
corresponds  to  an  harmonic  oscillator  degree  of  free¬ 
dom  \e ).  Then,  we  consider  the  composite  state 

l^>5  ®  \e)0, 

where  |^)s  =  ^ j  I <l>j)s  ls  the  state  of  the  system 

we  want  to  probe  and  \e)t  is  the  state  of  the  harmonic 
oscillator  in  the  coordinate  (jc-)  representation.  The 
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corresponding  state  in  the  momentum  ( p -)  represen¬ 
tation  is  denoted  \e)t. 

Assume  the  observable  (/-independent  Hermitian 
operator)  we  want  to  measure  is  A.  Then,  we  are 
interested  in  determining  s{^/\A\[l/)s  in  an  efficient 
way.  Suppose  that  we  know  how  to  implement  the 
unitary  operation  Us(t)  =  e-1^.  Following  Kitaev  we 
want  to  implement  the  following  conditional  evolution 

U  =  Y,\e)tAe\Us(t). 

t 

From  the  spectral  theorem  we  can  write  A  = 
T.j  Aj\4>j)ssi<Pj\-^n, 

M\(pj)s  ®  |0)o  =  ^Pu\<pj)s  <8>  \e), 

t 

=  4>j)s  ®  I e)t 

t 

=  I (Pj)s  ®  \Aj)t, 

where  |0)o  is  a  state  with  ( p  =  0)  zero  momentum. 
Basically,  the  conditional  evolution  U  is  a  momentum 
translation  operator  for  the  harmonic  oscillator  extra 
state.  Finally, 

^>s®|6)o  =  X>;l^>s®I^A- 

j 

Although  the  second  measurement  method  is  con¬ 
ceptually  simpler,  it  requires  approximately  imple¬ 
menting  the  ancillary  harmonic  oscillator,  the  condi¬ 
tional  evolutions  for  many  different  choices  of  r,  and  a 
more  complex  analysis  of  the  measurement  statistics. 
The  conditional  evolutions  can  be  simplified  some¬ 
what,  and  in  special  cases  (such  as  as  a  subroutine  of 
factoring)  become  very  efficient — see  [35]. 

4.4.  Measurement  noise  control 

The  quantum  physics  simulation  algorithm  de¬ 
scribed  above  is  approximate  and  the  output  is  noisy. 
In  order  to  properly  use  it,  we  need  to  have  explicit 
estimates  of  the  error  e  in  the  inferred  expectations 
given  the  noise  in  the  implementation.  Furthermore, 
the  effort  required  to  make  e  small  must  scale  polyno¬ 
mial^  with  1  /e.  There  are  three  sources  of  error  that 
need  to  be  considered.  The  first  is  associated  with  in¬ 
trinsic  noise  in  the  implementation  of  the  gates  due  to 
imperfections  and  unwanted  interactions.  The  second 


comes  from  the  discretization  of  the  evolution  opera¬ 
tor  and  the  use  of  the  Trotter  decomposition.  The  third 
is  due  to  the  statistics  in  measuring  the  desired  corre¬ 
lation  function  using  the  technique  given  above. 

4.4.1.  Gate  imperfections 

The  problem  of  gate  imperfections  can  be  dealt 
with  by  using  quantum  error  correction  [38,39]  and 
fault  tolerant  quantum  computation  [40-44].  Accord¬ 
ing  to  the  accuracy  threshold  theorem,  provided  the 
physical  gates  have  sufficiently  low  error,  it  is  pos¬ 
sible  to  quantum  compute  arbitrarily  accurately.  The 
fault  tolerant  computation  implements  unitary  opera¬ 
tions  and  measurements  on  encoded  qubits  with  over¬ 
heads  bounded  by  G( log* (I/O)  for  some  k.  This  ex¬ 
ponentially  efficient  convergence  implies  that  the  ef¬ 
fects  of  physical  noise  can  in  principle  be  ignored. 

4.4.2.  Discretization  error 

A  second  type  of  error  is  the  one  introduced 
by  the  discretization  of  the  evolution  operator.  This 
discretization  is  very  similar  to  the  one  used  in 
classical  simulation  of  dynamical  quantum  systems. 
It  is  possible  to  estimate  the  size  of  this  error  by  a 
detailed  analysis  of  the  discretization.  For  example, 
using  the  Trotter  approximation 

e-i(A/,+//2)Af  _  e-i//,Af/2e-i//2A/e— i//,A//2 

+  O((A03). 

The  coefficient  of  (A/)3  ~  —  i(//j  +  Hi)^ /6  can  be 
bounded  by  estimating  the  largest  eigenvalue  of  H\ 
and  Hj. 

4. 4. 3.  Measurement  statistics 

Our  technique  for  measuring  the  correlation  func¬ 
tion  (A(/)#(0))  requires  measuring  the  expectations 
of  unitary  operators  uJVy  associated  with  the  im¬ 
plemented  evolution.  In  most  cases,  the  operators  A 
and  B  are  a  sum  of  O^iua^b)  products  of  Pauli  matri¬ 
ces,  so  that  0(m/[)  U J’s  and  O(ms)  Vy's  arc  needed. 
This  means  that  the  expectation  is  a  sum  of  OfnAnis) 
random  variables  r;y,  where  \rjy\  ^  1.  To  assure 
that  the  statistical  noise  (given  by  the  standard  devi¬ 
ation)  is  less  than  c  it  suffices  to  measure  each  rjy 
0(mAtnB/c2)  times. 
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5.  Resonant  impurity  scattering 

5. 1.  Formulation  of  the  physical  problem 


Our  toy  problem  is  a  ring  of  n  equally-spaced 
lattice  sites  on  which  spinless  fermions  hop  to  nearest 
neighbor  sites  or  hop  onto  or  from  an  “impurity”  state. 
The  length  of  the  ring  is  L  =  na,  where  a  is  the 
distance  between  sites.  The  system  is  described  by  the 
Hamiltonian  (in  second  quantized  form) 


H  =  —T  ^(cjc,+i  +  cf+la)  +  (b^b 

i  =  l 

+  ^I2(c!b+btc0’ 


(22) 


i=l 


where  T  is  the  hopping  matrix  element,  e  is  the 
energy  of  the  localized  (impurity)  state,  and  V  is  a 
hybridization  energy.  As  usual,  b' s  and  c’s  are  fermion 
(anticommuting)  operators.  The  index  i  labels  the 
lattice  sites  (/?,  —  ia  is  the  lattice  site  position)  and 
strict  periodic  boundary  conditions  are  assumed,  i.e. 


(23) 


cf  =cf 
c/+n  ci  ‘ 

We  now  imagine  that  the  system  is  initially  pre¬ 
pared  in  the  zero  temperature  ground  state  of  the  ring 
in  the  absence  of  the  impurity.  Then,  at  time  t  =  0,  a 
fermion  is  injected  into  the  impurity  state.  After  the 
system  has  evolved  for  some  time  t ,  we  want  to  com¬ 
pute  the  probability  amplitude  that  the  evolved  state  is 
still  in  the  initial  state.  The  relevant  quantity  to  com¬ 
pute  is  (h  =  1  and  t  ^  0) 


G(O  =  (^(0)|fc(O^(0)|^(0)), 
b(t)  =  c'H'b(0)e~'Hl , 


(24) 

(25) 


where  the  initial  state  is  the  Fermi  sea  of  Ne  ^  n 
fermions 

Ne-  1 

|*(0))=|FS)=  n  cj.|0>.  (26) 

i=0 

|0)  is  the  vacuum  of  fermions  and 

4  =  ^Ee*'ieI  <27> 

The  wave  number  kj  is  determined  from  the  periodic 
boundary  conditions,  =  cj,  which  implies 


*7 


2izn  j 
L  ’ 


with  tij  an  integer. 


(28) 


There  is  no  unique  way  to  choose  the  set  of  tij ’s.  The 
common  convention  is  to  define  the  first  Brillouin  zone 
as 

71  ,  71 

- <k^-y  (29) 

a  a 


with  k  values  uniformly  distributed  in  this  interval 
with  spacing  A k  =  2tt/L. 


5.2.  Quantum  algorithm 


We  want  to  write  a  quantum  algorithm  that  allows 
one  to  compute  G(t).  To  this  end,  we  start  by  repre¬ 
senting  fermion  operators  in  terms  of  Pauli  matrices. 
Because  of  the  form  of  the  hybridization  term  a  most 
convenient  representation  is  the  following 


o 


2 

+  ’ 


from  which  the  following  mapping  results 

tfb  =  5(1  +  CT.1), 

4c*<  =  3(1 +^+2). 

clob+b^cko  =  i  (al  a2x  +  a'y  al ) . 
Therefore,  the  Hamiltonian  operator  reads 


/+ 2 


n- 1  1  n- 1 

2 H  =  €  -f  14-  eo]  -f  Ye^o? 

i=  0  J  1=0 

4-  V (crj o*  4- cty cty),  (30) 

where  £*  =  —2 T cos ka.  An  additional  simplification 
can  be  introduced  when  one  realizes  that  the  structure 
of  the  observable  to  be  measured  is  such  that 

b(t)  =  clHtb(  0)e~lHt  =  e,//rcrie_l///, 


where  H  is  given  by 


—  e 


£k0 


(31) 


(32) 


H  =  -0}  4-  ~f~Gz  +  Y {GX °X  +  CTyCTy), 

and,  therefore,  the  “string”  one  has  to  simulate  has 
length  equal  to  two  (it  involves  only  qubits  1  and  2) 


A(t)  =  b{t)b f  (0)  =  e^'crle-^'aj. 


(33) 
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If  we  were  to  transform  H  =  U  Hp\U t  unitarily  with 

U  =  UU  Q'HJpitj  an<^  n  a  ^n^te  integer  (UU^  =  1)  in 
such  a  way  that  both  Hp  i  and  Hpi  are  physical  Hamil¬ 
tonians,  then  the  simulation  would  be  straightforward. 
(We  call  this  type  of  mapping  a  physical  unitary  map¬ 
ping.)  For  our  two  qubit  case,  one  can  always  perform 
a  physical  unitary  mapping  with 


jj  gif  Gx  Q— i?' °y 


zuz  e'4w>  e' 4 
•  7T  _2  1_2  ;  7T  _  I  :  7T  _ 2 

^  g  '  ^  Oy  gl  J  (7,  Gy  g  Gj  q  1  4  Gy 


,  a  _  1  :  7T  _  2 


(34) 


///>!  =  ^(£-VA2  +  V2)^1 

+  ^(£  +  v/a2  + V2)ct2,  (35) 

with  £  =  (e  +  £*0)/2,  A  =  (€  —  £*o)/ 2,  and  cos#  = 
1  /Vi  +  «$2  with  5  =  (A  +  VA2  +  V2 )/  V. 

In  general,  such  a  constrained  transformation  is  not 
easily  realized  and  one  performs  a  Trotter  decomposi¬ 
tion 

=  [e*"*]"  =  +  Ofc2)]"  (36) 

where  //  =  //.  +  //(V  with  f/lv  =  y  (crjcr2  +  o^cr2) 
and  time  slice  s  =  t/M.  On  the  other  hand,  one  can 
easily  perform  a  physical  unitary  mapping  for  e'H*ys 


eiHxys  =  Ue{HplSu\  (37) 

where  Hp\  =  j (crj  -  crv2)  and 

=  ei?ffVi*<7'e-i?,T<<  (38) 


Finally,  the  “string”  one  has  to  simulate  with  the 
quantum  computer  is 

_  _  _  (39) 

S(s)  =  eiH'sUeiHf,‘sUf, 

and  G(t)  =  (A{i)). 


and  showed  how  this  problem  is  avoided  in  quantum 
computing  simulation.  The  evolution  of  quantum  com¬ 
puters  are  intrinsically  quantum  mechanical  and  this  is 
the  main  difference  with  a  classical  computer  that  al¬ 
lows  one  to  solve  the  sign  problem.  We  studied  sources 
of  errors  in  a  quantum  computer,  such  as  gate  imper¬ 
fections  and  the  expansion  of  the  evolution  operator, 
and  argued  that  they  would  not  open  a  back  door  to  a 
problem  similar  to  the  sign  problem. 

We  gave  a  very  general  definition  of  what  a  model 
of  quantum  computation  is.  In  particular  and  because 
of  our  particular  interest,  i.e.  the  simulation  of  fermion 
systems,  we  described  the  standard  and  the  fermionic 
models  (“Grassmann  Chip”).  These  arc,  of  course, 
not  the  only  ones.  Isomorphisms  of  *-algcbras  allow 
one  to  introduce  more  “esoteric”  models  [9].  Indeed, 
there  is  nothing  special  about  the  spinless  fermionic 
model  of  quantum  computation.  One  could  have  used 
a  “hard-core  boson”  model  which  admits,  in  principle, 
a  realization  in  terms  of  He4  atoms.  The  key  point  is 
the  implementation  of  the  physical  gates. 

Our  effort  focused  on  the  simulation  of  the  dy¬ 
namics  of  fermionic  quantum  systems.  However  other 
problems  can  be  of  interest:  the  thermodynamic  or 
ground  state  properties  of  a  Hamiltonian.  Even  if  one 
had  a  quantum  computer,  it  is  not  clear  how  to  use  it 
to  efficiently  compute  these  quantities.  On  the  other 
hand,  at  present,  no  proof  exists  showing  that  this  is 
not  possible. 

An  approach  that  in  principle  could  be  used  to 
compute  the  spectrum  of  a  Hamiltonian  H  (e.g.,  the 
ground  state)  or  expectation  values  of  arbitrary  observ¬ 
ables  is  the  adiabatic  “switching  on”  in  conjunction 
with  the  Gell-Mann-Low  theorem  [45]  of  quantum 
field  theory.  The  idea  simply  consists  of  introducing 
a  fictitious  Hamiltonian 


6.  Concluding  remarks 

We  investigated  the  implementation  of  algorithms 
for  the  simulation  of  fermionic  quantum  systems,  and 
gave  an  explicit  mapping  that  relates  the  usual  qubit 
of  a  quantum  computer  to  the  fermionic  modes  that 
we  want  to  simulate.  Our  attention  focused  on  the  so- 
called  sign  problem.  It  is  a  problem  appearing  in  at¬ 
tempts  to  simulate  classically  the  dynamics  of  quan¬ 
tum  systems.  We  reviewed  the  origin  of  this  problem 


H€(t)  =  H0  +  f€(t)Hu  (40) 

where  both  Ho  and  H\  are  time  independent  operators 
(H  =  Ho  4-  H\)  and  the  scalar  function  fe(t)  is  such 
that  limr_>±oo  /e(0  =  0  and  lim,_>o/e(0  =  1,  for 
an  arbitrary  adiabatic  parameter  e.  In  other  words, 
H€(t  =  0)  =  H  and  He(t  =  ±oo)  =  H0.  H0  is 
typically  an  operator  whose  spectrum  is  known,  e.g., 
an  arbitrary  bilinear  operator  representing  a  mean- 
field  solution  of  H  and  whose  eigenstates  can  be 
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straightforwardly  prepared  (let’s  call  it  |*o))-  The 
Gell-Mann-Low  theorem  asserts  that 

r  t/6(0, -oo)|*o)  l^o) 

lim - = -  (41) 

€-0  (*o|t/6(Of-oo)|*0>  <^ol^o) 

if  the  state  whose  limit  one  is  performing  admits  a  se¬ 
ries  expansion  in  a  coupling  parameter  characterizing 
the  strength  of  H\.  This  formal  device  generates  the 
eigenstate  adiabatically  connected  to  |<Z>o)-  The  the¬ 
orem  does  not  guarantee  that  if  |<Z>o)  is  the  ground 
state  of  Ho  then  Itf'o)  is  the  ground  state  of  H .  If  the 
conditions  of  the  theorem  are  satisfied  then  computa¬ 
tion  of  the  spectrum  of  H  is  straightforward.  To  our 
knowledge  this  approach  has  never  been  implemented 
in  practice. 

The  work  presented  here  is  only  a  first  step  in  a  pro¬ 
gram  investigating  the  simulation  of  quantum  systems 
using  quantum  computers.  We  have  given  a  rather  ex¬ 
plicit  algorithm  for  a  simple  problem  and  we  will  in¬ 
crease  the  complexity  of  the  problems  in  the  work  to 
come.  An  interesting  problem  would  be  to  provide  al¬ 
gorithms  to  test  for  superconductivity  in  systems  such 
as  the  Hubbard  model.  Such  simulations  using  classi¬ 
cal  computers  cannot  unequivocally  answer  this  im¬ 
portant  question  because  of  the  sign  problem,  but  a 
quantum  computer  could. 
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Abstract 

We  review  the  basic  ideas  behind  the  quantum  lattice  Boltzmann  equation  (LBE),  and  present  a  few  thoughts  on  the  possible 
use  of  such  an  equation  for  simulating  quantum  many-body  problems  on  both  (parallel)  electronic  and  quantum  computers. 
©  2002  Published  by  Elsevier  Science  B.V. 
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1.  Quantum  mechanics  and  fluids 

Intriguing  analogies  between  quantum  mechanics 
and  fluid  mechanics  have  been  pointed  out  since  the 
earliest  days  of  quantum  theory  [1].  The  orthodox 
tenet  is  that  these  analogies  are  purely  formal  in  char¬ 
acter  and  do  not  bear  upon  the  basic  physics  of  quan¬ 
tum  phenomena.  A  less-orthodox,  albeit  not  minor, 
stream  of  thought  insists  instead  that  quantum  me¬ 
chanics,  and  notably  Heisenberg’s  uncertainty  prin¬ 
ciple,  are  nothing  but  a  mirror  of  our  ignorance  of 
the  underlying  (hidden)  microscopic  physical  level. 
This  leads  to  the  puzzling  theory  of  ‘hidden  variables’ 
which  traces  back  to  Einstein  and  subsequently  to 
Bohm  and  others  [2].  It  is  not  our  intent  here  to  en¬ 
ter  this  fascinating  and  still  open  subject  [3].  We  turn 
to  a  practical  question  instead:  what  can  the  analogy 
do  for  us  in  terms  of  numerical  modeling  of  evolution¬ 
ary  quantum  mechanical  phenomena?  The  question  is 
legitimate  because,  regardless  of  its  philosophical  im¬ 
plications,  the  fluid  analogy  certainly  provides  an  in- 
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tuitive  and  physical  sound  basis  to  develop  numerical 
methods  for  time-dependent  quantum  mechanics.  In 
particular,  it  is  reasonable  to  ask  whether  the  advan¬ 
tages  brought  about  by  lattice  kinetic  methods  in  fluid 
dynamics  can — by  means  of  the  fluid  analogy — be  ex¬ 
ported  to  the  context  of  quantum  mechanics.  Before 
we  put  forward  our  discrete  kinetic  theory  version  of 
the  analogy,  it  is  useful  to  provide  a  cursory  survey  of 
the  main  ideas  behind  the  analogy  itself.  To  this  end,  a 
short  recap  of  basic  notions  of  quantum  mechanics  is 
in  order. 


2.  The  fluid  formulation  of  the  Schrodinger 
equation 

Let  us  begin  with  the  Schrodinger  equation  for 
a  non-relativistic  quantum  particle  of  mass  m  in  an 
external  potential  V(x): 

r  h2 

i hdtv  = - A  +  V(jc)  V,  (1) 

2  m 

where  ^(Jc,  t)  is  the  wavefunction  of  the  material  par¬ 
ticle.  Upon  multiplying  (1)  by  the  complex  conjugate 
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*/'*,  and  the  complex  conjugate  of  (1)  by  V  and  then 
subtracting,  we  obtain  the  following  set  of  fluid  equa¬ 
tions: 


d/P  +  daJa  =  0,  (2) 

dtJa  +  daP=  0,  (3) 

where,  by  using  the  eikonal  representation  ^  = 


pl/2eW; 

P=  I'J'I2. 

(4) 

h 

Ja  =  =  pua, 

m 

(5) 

fu2  V 

P  —  P\-T-  - 1 - ), 

\  2  m  m  ) 

(6) 

and 

vG__*2  v/2 

2m  p1/2 

(7) 

is  the  famous  quantum  potential  advocated  by  Bohm 
and  coworkers  to  support  the  picture  of  quantum  me¬ 
chanics  as  an  intrinsically  non-local  description  of 
the  microscopic  world  [4].  This  configures  quantum 
matter  as  an  ideal  ( inviscid ,  dissipationless),  irrota- 
tional  compressible  fluid.  The  inviscid  character  of  the 
quantum  fluid  stems  from  the  reversible  nature  of  the 
Schrodinger  equation,  a  diffusion  equation  in  imagi¬ 
nary  time.  So  much  for  the  analogy  in  the  continuum. 
What  about  the  discrete  lattice  world?  Interestingly 
enough,  this  analogy  becomes  even  more  compelling 
once  transposed  into  the  language  of  the  lattice  world. 
In  fact,  the  lattice  formulation  naturally  calls  for  an  an 
“upgrade”  from  the  non-relativistic  Schrodinger  equa¬ 
tion  to  its  relativistic  associate,  the  Dirac  equation. 
Symbolically,  the  analogy  goes  as  follows  ( DE :  Dirac 
equation,  SE:  Schrodinger  equation,  LBE:  Lattice 
Boltzmann  equation,  NSE:  Navier-Stokes  equation): 

DE  -*  SE ,  (8) 

LBE  NSE.  (9) 

2.  ].  Fluid  formulation  of  relativistic  quantum 
mechanics 

To  unfold  this  analogy,  it  proves  expedient  to  cast 
the  Dirac  equation  into  a  form  where  all  streaming 
matrices,  known  as  Weil  matrices,  become  real.  This 
is  the  so-called  Maiorana  form :  In  a  compact  four¬ 
dimensional  notation,  this  reads 


[w?kdv]xl'k  =  iMjk\l/k,  n  =  0,3  (10) 

with 

W%=SJk, 

w}k=Pji"  K  =  -aW 

M  jk  =  ~mayJk  +  qVSjk  +  Aajk  Ja , 

where  all  matrices  have  the  standard  meaning.  Here 
JflAtl=qV  +  JaAa  is  the  interaction  of  the  elemen¬ 
tary  charge  q  with  an  external  electromagnetic  field 
described  by  the  4-vector  potential  ( V ,  Aa). 

A  scalar  product  of  Eq.  (10)  with  \//*  yields  the 
desired  set  of  continuity  equations: 

atPj  +  9aJj=Sj,  7  =  1,4  (11) 

where  pj  —  is  the  partial  density  of  the  7  th  fluid, 

Jj  =  the  corresponding  current  density,  and 

Sj  =  \\//* M jk^k  is  a  “chemical”  source  term  transfer¬ 
ring  mass  across  the  different  components  of  the  rel¬ 
ativistic  mixture.  Note  that  in  the  above  expressions 
only  the  index  k  is  summed  upon.  Unitarity,  read  norm 
conservation,  implies  Yj  =  Yljk  ^j^jk^k  =  0. 
This  is  automatically  secured  by  the  antihermitian 
character  of  the  mass  matrix:  My  +  M*k  =  0.  As 
promised,  the  fluid  analogy  comes  by  more  naturally 
than  in  the  non-relativistic  case,  because  the  Dirac 
equation  only  involves  first  order  derivatives.  Another 
pleasing  feature  is  that  the  external  interaction  is  easily 
accommodated  into  a  formal  redefinition  of  the  mass 
matrix,  without  compromising  the  local  nature  of  the 
theory.  The  fluid  interpretation  of  the  Dirac  equation 
is  equally  transparent:  four  types  of  spinning  particles 
stream  in  space  and,  once  on  the  same  space-time  lo¬ 
cation,  they  interact  via  the  “scattering  matrix” 

A  qualitative  difference  with  classical  particle  motion 
is  apparent,  though.  In  a  classical  fluid,  particles  do  not 
“mix”  during  the  streaming  phase.  A  type-1  particle  at 
location  jc  at  time  t  with  speed  v  propagates  to  x  -f  v  dt 
at  time  t  +  dr  and  it  is  still  entirely  of  type  1 . 

A  relativistic  particle  however  undergoes  mixing 
during  free  propagation,  because  its  spinning  mo¬ 
tion  implies  a  rotation  around  the  direction  of  motion 
which  mixes  up  the  four  spinorial  components.  This  is 
why  the  streaming  matrix  is  generally  non-diagonal, 
echoing  the  fact  that  spin  is  not  an  ordinary  vector. 
This  suggests  that  the  discrete  space-time  of  a  rela¬ 
tivistic  particle  should  be  represented  by  a  ‘hypemet- 
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ted  lattice’  in  which  each  link  is  made  up  of  four  dis¬ 
tinct  but  communicating  channels,  one  per  spinorial 
state.  This  “Hypemetted  Lattice  Theory”  is  less  of  a 
joke  than  it  seems.  It  has  been  recently  realized  that 
lattice  formulations  of  field  theory  based  upon  spin¬ 
ning  particle  motion  may  offer  potential  advantages 
over  more  popular  techniques  such  as  path-integration 
[5].  This  is  because  in  quantum  lattice  models  “ instead 
of  seeking  discretized  versions  of  the  Hamiltonian  or 
the  Lagrangian,  a  discretized  version  of  the  evolution 
operator  is  introduced ”  [5].  In  fact,  what  this  author 
finds  is  that  “the  rotation  group,  the  Lorentz  group  and 
spin  emerge  automatically  in  the  continuum  limit  from 
unitary  dynamics  on  a  cubic  lattice ”.  The  reader  fond 
of  more  details  is  directed  to  the  original  reference. 

2.2.  Dirac  to  Schrodinger:  the  adiabatic 
approximation 

As  noted  in  [6],  the  way  the  Schrodinger  equation 
is  obtained  as  a  long  wavelength  (low  energy)  limit 
of  the  Dirac  equation  involves  a  sort  of  adiabatic 
approximation  which  is  formally  very  similar  to  the 
low-Knudsen  adiabatic  expansion  taking  the  (lattice) 
Boltzmann  equation  into  the  Navier-Stokes  equations. 
The  formal  parallel  emerging  from  this  analogy  is 

Kn  =  Ip/ Im  ^  ft  =  v/c,  (12) 

where  lfl  is  the  particle  mean  free  path,  Im  a  typical 
coherence  length  of  the  macroscopic  fluid  and  /  is 
the  relativistic  particle  to  light  speed  ratio.  Detailed 
calculations  can  be  found  in  the  original  reference  [6] 
and  need  not  be  repeated  here. 

The  relativistic  motion  implies  that  any  particle  of 
momentum  pa  is  invariably  associated  with  an  antipar¬ 
ticle  with  opposed  momentum  —  pa.  The  symmetric 
combination  of  these  two  gives  rise  to  a  smooth,  emer¬ 
gent  field,  whereas  the  antisymmetric  combina¬ 
tion  defines  a  low  amplitude,  high-frequency  mode 
which  decouples  from  the  system  dynamics  in  the 
limit  p  — ►  0.  The  scenario  is  exactly  the  same  as  the 
adiabatic  approximation  in  kinetic  theory,  with  a  key 
difference.  Kinetic  theory  describes  dissipative  phe¬ 
nomena  in  which  adiabatic  elimination  wipes  out  the 
initial  conditions,  the  transient  modes  die  out,  never 
to  return.  Quantum  mechanics  is  reversible,  and  fast 
modes  never  die  out:  they  just  oscillate  so  fast  that 


any  observation  on  timescales  longer  than  their  pe¬ 
riod  of  oscillation  simply  overlooks  them.  But  they  are 
still  there  and  more  resolved  (higher  energy)  measure¬ 
ments  could  always  bring  them  back  again.  Note  that 
it  is  the  fast  mode,  not  the  antiparticle  mode  that  fades 
away;  the  particle-antiparticle  twin-link  does  not  dis¬ 
solve  even  in  the  low  energy  limit. 

Another  interesting  remark  concerns  the  symmetry 
breaking  induced  by  a  non-zero  mass  m.  If  m  is 
made  zero  the  up  and  down  walkers  do  not  see  each 
other  and  go  across  with  no  interaction,  the  result 
being  the  wave  equation  for  photons.  Manifestly  this 
is  a  singular  limit  which  cannot  be  described  by 
the  Schrodinger  equation  (diffusion  coefficient  goes 
to  infinity).  Any  non  zero  mass  causes  “collisions” 
which  slow  down  the  wavepackets  and  confer  them  a 
subluminal  speed  v  <  c  as  it  befits  material  particles. 

2.3.  The  interacting  case 

Interactions  with  an  external  or  self-consistent 
fields  are  readily  included  by  a  minor  extension  of  the 
“collision  operator”.  They  read  as  follows: 

dtu\2  -dzu\,2=md2,\  +  igd2.i,  (13) 

dtd\,2  +  M  1,2  =  -mu2A  +  1.2,  (14) 

where  g  =  eV/h  is  the  coupling  frequency  of  the 
potential.  Self-consistent  potentials,  such  as  those 
arising  in  connection  with  the  non-linear  Schrodinger 
equation,  are  easily  accommodated  by  making  g  a 
function  of  the  local  density  u2  -f  d2. 

3.  The  quantum  lattice  Boltzmann  equation 

We  are  finally  in  the  position  to  reformulate  the 
basic  analogy  in  quantitative  terms.  This  is  based 
on  the  following  position:  The  4-spinor  (Jc ,  r)  = 
\l/(x,Sj,t)  is  identified  with  a  complex  discrete  par¬ 
ticle  distribution  f  (x ,  t)  =  /(jc,  5/,  t).  The  analogy  is 
tantalizing,  but  a  minute’s  thought  reveals  two  severe 
flaws: 

(1)  While  the  4-spinor  \//j  (we  consider  spin  1/2 
throughout)  has  always  four  components  in  any 
dimensions,  the  discrete  population  /•  is  a  set  of  b 
real  functions  with  b  a  sensitive  function  of  space 
dimensionality. 


320 


S.  Sued  /  Computer  Physics  Communications  146  (2002)  317-323 


(2)  While  LBE  streaming  is  always  diagonal  in  mo¬ 
mentum  space,  the  three  Weil  matrices  can  not  be 
simultaneously  diagonalized. 

Both  problems  are  intimately  related  to  the  quantum 
nature  of  the  spin  variable.  Fortunately,  there  is  a 
way  out.  As  observed  in  [6]  in  the  limit  of  ‘small’ 
timesteps,  actually  much  shorter  than  the  inverse 
Compton  frequency  a;"1,  both  flaws  can  be  circum¬ 
vented  by  decomposing  the  three-dimensional  particle 
motion  into  a  sequence  of  three  one-dimensional  mo¬ 
tions  along  the  coordinate  directions  x ,  y,  z.  The  tech¬ 
nical  key  to  achieve  this  task  is  a  well  known  tool-of- 
the-trade  in  computational  fluid  dynamics:  “Operator 
Splitting’’.  The  main  use  of  operator  splitting  in  Com¬ 
putational  Fluid  Dynamics  is  to  handle  3D  problems 
as  a  sequence  of  lower  dimensional  ones.  In  quantum 
field  theory,  a  very  similar  technique  goes  under  a  dif¬ 
ferent  name:  “Trotter  formula”:  qa  =  (QAln)n  with  n 
integer  and  A  any  ‘reasonable’  operator.  Consider  the 
formal  solution  to  the  Dirac  equation  for  a  massless 
particle  (the  collisional  operator  plays  no  role  at  this 
stage): 

*j(  X11  +  dr„)  =  [ed'E'3‘=°  Wi^‘\pk{xIL).  (15) 

Manifestly,  the  propagator  taking  the  wavefunc- 
tion  from  x =  (jc ,  y,z,r)  to  xtl  +  dxfl  =  (x  + 
d*,y  +  dy,z  +  dz,t  +  dr)  is  the  direct  product 

of  three  one-dimensional  partial  propagators  Pa  = 

dfia,+w?.aa] 

e  Jk  ,  a  =  x,  y,  z  (no  summation  upon  a  im¬ 
plied).  This  is  the  natural  consequence  of  the  additiv¬ 
ity  of  the  streaming  operator.  This  expression  is  a  good 
starting  point  for  “conventional”  numerical  treatment 
of  the  Dirac  equation  [7],  but  is  definitely  unsuitable  to 
a  quantum  LBE  formulation  because  spinorial  states 
get  mixed  during  the  propagation  step,  something  that 
would  not  occur  to  a  classical  particle. 

Therefore,  a  naive  application  of  operator  splitting 
is  not  viable. 

However,  we  can  argue  that  we  do  not  need  to  work 
with  the  same  representation  of  the  Dirac  equation 
during  the  three  separate  streaming  steps.  As  long  as 
we  are  able  to  develop  a  recipe  securing  uniqueness  of 
the  representation  in  x[L  and  x[L  -F  At,*,  we  are  free  of 
choosing  the  representation  that  better  fits  our  needs. 
The  idea  is  to  perform  each  ID  partial  streaming 
in  the  representation  where  the  corresponding  Weil 
matrix  is  diagonal.  In  practice,  one  propagates  along 


one  direction,  say  jc,  then  ‘rotates’  the  system  so 
as  to  diagonalize  the  Weil  matrix  along,  say,  y,  so 
that  propagation  along  y  can  be  performed  like  for 
a  classical  particle,  and  finally  ‘rotates  back’  the 
propagated  solution  at  (jc  +  dje,  y  -f  dy,  t  +  dr).  New 
errors  are  introduced  in  the  numerical  treatment,  but 
we  shall  argue  that  they  are  0(dr2),  namely  within  the 
general  accuracy  of  the  LBE  method. 

The  quantum  LBE  bears  many  similarities  with 
other  quantum  lattice  schemes  discussed  in  the  recent 
[5,8,9]  and  not  so  recent  [10]  literature.  What  sets  it 
apart  from  all  these  schemes  is  the  fact  of  insisting 
on  a  diagonal  representation  of  the  Weil  matrices, 
so  as  to  retain  the  notion  of  classical  trajectories 
as  much  as  we  can.  In  fact,  the  “turn”  operator  R 
can  formally  be  interpreted  as  an  “internal  scattering” 
between  particle-antiparticle  states  [9],  thus  leaving 
the  concept  of  quantum  trajectory  still  well  defined, 
although  in  a  generalized  sense.  In  a  pictorial  sense 
[10],  we  might  say  that  while  classical  particles 
just  “Stream  and  Collide”,  quantum  particles,  like 
swimmers,  need  a  somersault  before  they  can  turn  in 
space:  they  “Stream,  Turn  and  Collide”!  The  ‘Turn’ 
step  is  a  necessity  induced  by  the  internal  structure  of 
the  relativistic  particle. 

Leaving  the  details  to  the  original  work,  here 
we  simply  report  the  final  result  for  a  pair  of  ‘up’ 
and  ‘down’  walkers  in  one-dimension.  Upon  using 
a  Cranck-Nicholson  time-marching  procedure  (secur¬ 
ing  unitarity  of  the  numerical  scheme),  the  quantum 
LBE  takes  the  following  form: 

u(z  4-  dz,  t  -f  d t)  =  Au(z,  t)  +  Bd(z ,  f),  (16) 

d(z  -  dz,  t  4-  dr)  =  Ad(z,  t)  -  Bu(z,  0,  (17) 

where 

A  —  '-W  . 

1  -  £2/4  -  i g 

B  — - — - , 

l-fl/4-i  g 

n=m2-g2. 

A  few  comments  are  in  order. 

First,  with  g  =  0  (no-interaction),  implicit  time 
marching  translates  into  a  mere  redefinition  of  the 
particle  mass  m  — ►  m!  —  m/{\  —  m2/ 4).  By  rein¬ 
stating  the  time-step  A r,  it  is  easily  recognized  that 


(18) 

(19) 

(20) 
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m'  ->  m  in  the  limit  At  ->  0,  which  means  that  quan¬ 
tum  LBE  fulfills  the  requirement  of  numerical  con¬ 
sistency.  Large  timesteps  mAt  >  1  lead  to  unphys¬ 
ical  results,  as  it  is  to  be  expected  since  the  natural 
Compton  frequency  m  (in  atomic  lattice  units  h  =c  = 
At  =  Ax  =  1),  is  no  longer  resolved.  Simple  alge¬ 
bra  also  shows  that  quantum  LBE  is  unconditionally 
stable  and  norm-preserving  (the  all-important  unitar- 
ity  condition).  This  is  fairly  remarkable  for  an  explicit 
numerical  scheme  [11],  and  ultimately  traces  back  to 
the  (implicit)  lightcone  discretization  hidden  behind 
the  quantum  LBE,  Eq.  (16).  Finally,  note  that  at  no 
point  in  our  treatment  did  we  need  to  care  about  strin¬ 
gent  symmetry  requirements:  apparently  a  simple  cu¬ 
bic  lattice  is  good  enough  to  our  purpose.  This  prob¬ 
ably  relates  to  the  diagonal  nature  of  the  quantum- 
mechanical  pressure  tensor  and  to  the  fact  that,  unlike 
fluid  dynamics,  the  theory  is  not  self-interacting.  Fi¬ 
nally,  we  observe  that  quantum  LBE  is  as  computa¬ 
tionally  lean  and  amenable  to  parallel  processing  as  an 
explicit  scheme  can  be. 

All  in  all,  a  good  set  of  credentials  for  a  numerical 
scheme. 

4.  Numerical  tests 

The  quantum  LBE  scheme  has  been  validated  on 
a  series  of  one-dimensional  textbook  calculations, 
including 

(i)  free  particle  propagation, 

(ii)  harmonic  oscillator, 

(iii)  scattering  from  a  rectangular  barrier  [12]. 

In  addition,  the  scheme  has  also  been  demonstrated 
for  simple  cases  of  non-linear  Schrodinger  equations 
of  direct  relevance  to  Bose-Einstein  condensation 
[13]  (as  an  example,  see  Fig.  1).  These  tests  pro¬ 
vide  evidence  of  the  viability  of  the  quantum  LBE  in 
one-dimension.  The  scheme  performs  efficiently  and, 
what’s  more,  provides  stability  and  unitarity  at  a  time, 
a  very  valuable  property  for  an  explicit  scheme.  As  we 
said,  this  is  related  to  the  peculiar  light-cone  space- 
time  marching  technique  inherent  to  quantum  LBE. 
Higher-dimensional  versions  akin  to  the  quantum  LBE 
discussed  here  have  been  developed  systematically  by 
Boghosian  and  coworkers  [8]. 


5.  The  quantum  A-body  problem 

In  this  section  we  shall  explore  the  question  of 
whether/what  the  lattice  techniques  discussed  so  far 
can  bring  any  new  insight  into  the  problem  of  solving 
the  Schrodinger  equation  for  a  collection  of,  say,  A 
particles  (quantum  A-body  problem): 

N 

\hd,<t>  =  ^[— A„  +  V  (X„)]<f>,  (21) 

>1=1 

where  Xn  =  (xn,yn,zn)  is  the  spatial  coordinate  of 
the  ;?th  particle,  d>(X\  . . .  Xjy)  the  A-body  wavefunc- 
tion  and  V  the  interparticle  potential,  typically  in  a 
two-body  format  V(Xn)  =  Zm>n  V(| Xn  -  Xm\).  It 
has  been  recently  pointed  out  [9]  that  quantum  lattice 
algorithms  constitute  excellent  candidates  as  numer¬ 
ical  schemes  for  quantum  computers.  In  the  A-body 
quantum  LBE,  each  quantum  particle  is  represented 
by  bG  walkers,  b  being  the  coordination  number  of 
the  lattice,  namely  the  number  of  discrete  momen¬ 
tum  states  attached  to  each  lattice  site.  These  walkers 
move  around  according  to  a  fictitious  microdynam¬ 
ics  whose  macroscopic  limit  is  precisely  the  A-body 
Schrodinger  equation. 

What  would  this  A-body  quantum  LBE  algorithm 
look  like? 

“Simply”  evolve  A  replicas  of  the  single-particle 
quantum  LBE  scheme  and  tie  them  up  together  via 
a  two-body  potential  collecting  the  sum  of  all  con¬ 
tributions  V"  =  Ylg'  Jim  v(xg  ~  x'p  at  each  8ivcn 
site  X".  If  one  does  not  insist  on  the  idea  of  a 
particle  generalized-trajectory,  and  turns  instead  to  a 
‘information-network’  picture,  a  generic  quantum  lat¬ 
tice  algorithm  would  take  the  form  of  a  first-order,  ex¬ 
plicit,  non-local,  map  for  the  complex  array  <Z>, : 

Vj(X,t)  =  YJTjMX-  VkAt,t-  AO,  (22) 

k 

where  V*  scans  the  3A-dimensional  neighborhood 
of  Xn  =  ( xn,yn,zn ),  n  =  1,  A  and  Tjk  is  the  com¬ 
plex  transfer  matrix  fulfilling  the  unitarity  condition 
TjiTik  =  Sjk .  The  kinetic  energy  operator  is  sweet 
since  any  walker  in  a  given  single-particle  state  can  be 
moved  independently  of  the  others,  resulting  in  a  lin¬ 
ear  0(bG  A)  complexity.  Unfortunately,  the  two-body 
long-range  potential  generates  a  daunting  quadratic 
complexity,  ( bGN )2,  to  say  nothing  of  the  ( bG)N  re¬ 
quirement  in  computer  storage  ....  The  scheme  meets 
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with  a  “exponential  complexity  wall”  which  rules  out 
any  possible  use  of  conventional  electronic  computers 
for  more  than  a  few  hundreds  particles  [17].  Although 
this  statement  can  probably  be  challenged  by  modern 
multiscale  techniques  (Achi  Brandt,  private  communi¬ 
cation),  we  shall  assume  that  such  exponential  com¬ 
plexity  is  indeed  beyond  electronic  computation  ca¬ 
pabilities.  This  brings  us  back  to  quantum  computers. 
Since  the  matter  of  solving  the  A-body  Schrodinger 
equation  in  full  on  a  quantum  computer  has  been  de¬ 
scribed  in  the  existing  literature,  here  we  shall  take 
a  different  path,  and  discuss  how  efficient  real-space 
single-particle  quantum  solvers  may  contribute  to  ad¬ 
vancing  the  A-body  frontier  without  solving  the  full 
N-body  Schrodinger  equation.  Incidentally  we  note 
that  this  is  of  actual  interest  not  only  for  current  elec¬ 
tronic  parallel  computers,  but  hopefully  also  for  actual 
software  emulators  of  quantum  computers  [14,15]. 

5. 1.  Quantum  LBE  and  Density  Functional  Theory 

As  previously  discussed,  numerical  algorithms  for 
the  quantum  many-body  wavefunction  are  very  hard 
(to  say  the  least)  on  electronic  computers.  Many  ways 
out  have  been  developed  to  cope  with  this  problem, 
including, 

(i)  Quantum  Monte  Carlo  techniques  [18], 

(ii)  Multiscale  methods  [19], 

(iii)  Effective  one-body  theories. 

In  this  paper,  we  shall  be  concerned  with  option 
(iii). 

Effective  one-body  theories  developed  in  the  last 
fourteen  years  permit  to  learn  a  great  deal  about  the 
properties  of  quantum  many-body  systems  without 
ever  invoking  the  use  of  many-body  wavefunctions. 
Particularly  successful  in  this  respect  is  the  famous 
Density  Functional  Theory  developed  in  the  60’s 
by  Hohenberg-Kohn  and  Kohn-Sham  [16,17].  The 
core  idea  of  Density  Functional  Theory  is  that  the 
ground  state  of  a  many-electron  wavefunction  (nuclei 
are  regarded  as  classical  particles  on  account  of 
their  higher  mass)  is  uniquely  determined  by  the 
electronic  density  n(x)  =  .  |0y|2(*),  where  are 

one-particle  orbitals.  The  ground-state  energy  can 
then  be  obtained  by  summing  up  the  single-particle 


orbital  energies  obtained  by  solving  the  Kohn-Sham 
equations: 

#ks  4>j  =  Ej<pj,  (23) 

where  the  Kohn-Sham  Hamiltonian  consists  of  four 
contributions 

h2 

HKS  =  “  —  A j  +  VextU) 

2m 

+  e2  ff/fdy+  Vexln].  (24) 

J  \y-x\ 

The  first  two  contributions  are  the  usual  kinetic  en¬ 
ergy  and  external  potential  operators,  the  third  one  re¬ 
lates  to  the  self-consistent  Hartree-Fock  potential.  Fi¬ 
nally,  the  fourth  one  is  an  effective  ‘exchange’  energy 
functional  which  collects  the  effects  of  A-body  inter¬ 
actions.  The  idea  is  that  an  effective  functional  of  the 
electron  density  exists  such  that  the  ground  state  en¬ 
ergy  of  a  fictitious  system  of  independent  electrons 
moving  in  such  a  potential  is  exactly  the  same  ground- 
state  energy  of  the  interacting  system!  Describing  how 
such  a  magic  comes  about  is  certainly  beyond  the 
scope  of  this  work.  Here,  we  shall  simply  remark 
that  Density  Functional  Theory  heavily  leans  on  the 
intuitive  picture  of  a  quantum  many-body  system  as 
a  backbone  of  ions  tied  up  together  by  a  very  mo¬ 
bile  electronic  fluid.  In  this  respect,  it  certainly  puts 
a  premium  on  efficient  real-space  solvers  for  the  one- 
particle  (non-linear)  Schrodinger  equation,  both  in  the 
time-independent  (ground-state)  and  time-dependent 
(excited  states)  form.  A  practical  scheme  which  could 
be  implemented  today  on  either  electronic  or  quantum 
computer  emulators  is  briefly  outlined  in  the  follow¬ 
ing. 

Consider  the  task  of  solving  a  set  of  A  effective 
time-dependent,  one-particle,  Kohn-Sham  equations 
coupled  via  an  effective  potential  VksIp]- 

i  hdt(pj  =  //ks0/ •  (25) 

Since  the  LBE  grid  is  uniform,  the  non-local  Hartree- 
Fock  potential  is  best  turned  into  the  correspond¬ 
ing  Poisson  problem  A  Vhf  =  n,  which  is  efficiently 
solved  by  standard  methods  such  as  rapid  ellip¬ 
tic  solvers  or  Fast-Fourier  techniques.  The  exchange 
functional  is  local  and  can  therefore  be  handled  by 
the  same  procedure  already  discussed  and  tested  for 
Bose-Einstein  condensation. 
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1  I -  I  I - > 

2  | -  |  | - > 

3  | -  |  | - > 

4  | -  |  | - > 

t  t+dt 

Sketch  1.  Parallel  solution  of  the  set  of  N  Kohn-Sham  equation 
(TV  =4).  The  double  line  =  indicates  the  serial  phase  in  which 
each  slave  processor  forwards  its  partial  density  to  the  master  and 
subsequently  receives  the  effective  potential  to  initiate  the  next  step. 


64  128  192  256  320  384  448  512 
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Fig.  1.  Probability  density  /)+(z)  at  t  =  0  and  t  =  25,000 
for  the  case  of  a  self-consistent  potential  of  the  form:  V(z)  = 
—  Zq)2  +  VnP-  The  flattening  of  the  wavefunction  due  to 
self-consistent  interactions  is  well  visible.  The  superscript  -I-  indi¬ 
cates  the  slow,  hydrodynamic  mode. 

The  system  of  Eqs.  (25)  is  also  particularly  well- 
suited  to  parallel  computing.  In  fact,  each  of  these 
equations  can  be  advanced  concurrently  in  time.  Upon 
completing  a  single  time-step  calculation,  each  proces¬ 
sor  forwards  its  partial  density  pj(t  +  df)  =  \(j)2\ (f  + 
df)  to  a  master  processor  whose  task  is  to  collect  all 


contributions,  form  the  effective  potential  Vfcs  and 
send  it  back  to  each  processor  to  initiate  the  next  time 
step.  This  process  can  be  performed  fairly  efficiently 
on  a  electronic  parallel  since  it  entails  a  very  lean 
communication-to-computation  ratio.  For  instance,  a 
parallel  computing  consisting  of  P  =  N  processors 
would  solve  the  N-body  problem  in  (slightly  more 
than)  the  same  time  it  takes  a  serial  one  to  solve  the 
single-body  quantum  equation.  Since  the  same  state¬ 
ment  applies  to  a  single-processor  quantum  computer, 
one  might  dream  of  quantum-computers  applications 
of  paramount  scientific  problems,  such  as  electronic 
structure  calculations  of  large  molecules  of  biological 
interest. 
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Abstract 

We  demonstrate  a  strategy  for  implementation  a  quantum  full  adder  in  a  spin  chain  quantum  computer.  As  an  example,  we 
simulate  a  quantum  full  adder  in  a  chain  containing  201  spins.  Our  simulations  also  demonstrate  how  one  can  minimize  errors 
generated  by  non-resonant  effects.  ©  2002  Published  by  Elsevier  Science  B.V. 

PACS:  03. 67. Lx;  03.67.-a;  76.60.-k 
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1.  Introduction 

A  full  adder  is  a  basic  component  of  a  conventional 
computer  and  a  welcome  asset  for  quantum  comput¬ 
ers.  In  particular,  the  Shor  quantum  algorithm  requires 
modular  exponentiation,  f(x)  =  ax  (mod  V),  which 
cannot  be  computed  without  a  quantum  full  adder.  The 
question  arises:  What  is  a  full  adder  in  a  quantum  com¬ 
puter?  A  quantum  computer  operates  on  a  superposi¬ 
tion  of  numbers  simultaneously.  What  is  not  clear  is 
what  it  means  to  add  two  superpositions.  One  defin¬ 
ition  of  a  full  adder  in  a  quantum  computer  is  that  a 
full  adder  is  a  gate  which  adds  a  given  number  to  a 
superposition  of  numbers.  This  full  adder  must  simul¬ 
taneously  add  a  definite  number,  “a”,  to  all  numbers, 
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which  are  in  a  superposition  in  a  quantum  com¬ 
puter  register.  The  basic  idea  for  the  full  adder  is  well 
known  [1,2].  However  it  is  not  clear  how  to  implement 
this  idea  in  a  many-qubit  spin  quantum  computer.  It 
is  even  less  clear  how  to  simulate  accurately  a  quan¬ 
tum  full  adder  having  a  large  number  of  qubits  using  a 
conventional  computer.  One  must  understand  the  role 
of  non-resonant  effects  and  how  to  minimize  them. 
One  must  also  know  the  structure  of  the  error  states 
caused  by  non-resonant  effects.  Our  paper  provides 
an  answer  for  these  issues.  To  implement  a  quantum 
full  adder,  we  propose  to  use  a  sequence  of  electro¬ 
magnetic  7r -pulses  on  a  spin-chain  quantum  computer. 
We  have  simulated  the  dynamics  of  this  quantum  full 
adder  with  201  spins,  on  a  conventional  computer;  an¬ 
alyzed  unwanted  non-resonant  effects;  and  determined 
the  structure  of  the  error  states  and  ways  to  reduce 
non-resonant  effects. 
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Table  1 

Table  of  binary  addition 


a 

b 

c 

5' 

cf 

ah  ©  ac  ©  he 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

1 

0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

0 

1 

1 

0 

1 

1 

1 

0 

1 

0 

1 

1 

1 

1 

1 

1 

0 

1 

0 

1 

1 

1 

1 

1 

2.  A  classical  full  adder 

A  classical  full  adder  operates  with  an  input  of 
two  addend  bits,  “a”  and  “b”,  and  a  carry  bit,  “c”. 
(See  Table  1.)  In  Table  1,  s'  and  c'  are  the  output 
sum  and  the  carry-over,  respectively.  The  sum,  s',  can 
be  easily  expressed  as  s'  =  a  ©  b  ©  c  (where  0  is 
addition  modulo  2).  The  formula  for  the  carry-over 
is  more  complicated.  One  can  see  from  Table  1  that 
c'  =  ab  0  ac  0  be. 


3.  A  quantum  full  adder 

We  describe  the  basic  idea  of  a  quantum  full  adder 
first  suggested  in  [1].  The  full  adder  quantum  gate 
(F)  depends  on  the  value  of  the  bit  “a”.  If  a  =  0,  the 
quantum  computer  applies  the  gate  F(0).  If  a  =  1,  it 
applies  F(  1).  The  F-gate  can  be  expressed  in  terms  of 
the  Control-Not  (C,*)  and  Control-Control-Not(C/jtp) 
gates.  A  Cj k -gate  changes  the  value  of  the  target  qubit, 
“k”,  if  the  control  qubit,  “i”,  has  the  value  “1”.  A  CucP- 
gate  changes  the  state  of  the  target  qubit,  “p”,  if  both 
qubits,  “i”  and  “k”,  have  the  value  “1”, 

Ci*|... =  |  ...«i  ...(«/  ©«*)*••■). 

Cikp\ •  ■  •  «/ . .  .n/c .  ..ttp . . .)  (1) 

=  |...n; . .  .«*  . . .  [(«/«*)  ©  np]p . . 

We  shall  also  use  the  Not-gate  which  can  be 
designated  by  C/.  It  changes  the  value  of  the  ith  qubit 
independent  of  the  values  of  all  other  qubits.  Using 
these  gates,  we  can  transform  the  state,  |03C2^i ),  into 
the  state,  W^bx), 

F}2l(a)\0iC2b\)  -  \c'3s'2bi),  (2a) 


where  “a”  and  “b”  are  the  addend  bits  and  s'  is  their 
sum,  c  and  c'  are  the  input  and  output  carry  bits. 
Namely, 

^321(0)  =  C12C123, 

(2b) 

/r32l(l)  =  C12C123C2C23. 

We  use  the  convention  that  the  right  gate  acts  first. 
Let  us  check  the  action  of  the  full  adder,  F321,  for 
example,  for  a  =  b  =  c  =  1 .  We  use  F321  (1): 

C23|03C2M=C23|03l2ll>  =  |l3l2ll>.  (3) 

In  (3),  the  second  control  qubit  has  the  value  “1”.  So, 
the  third  target  qubit  changes  its  value  from  “0”  to  “  1 
Next, 

C2 1 1 3 1 2 1 1  >  —  I I3O2  1 1  >» 

C123II3O2I1)  =  I I3O2 1 1 ),  (4) 

C12I I3O2 1 1  >  =  1 13 12 1 1)  =  |c352^l)‘ 

Thus,  in  accordance  with  the  Table  1 ,  we  obtain, 
b  =  1,  s'  =  1,  c  =  1. 

4.  A  spin  chain  quantum  computer 

An  atomic  chain  quantum  computer  based  on 
triplets,  ABC  ABC  ABC  . . .,  was  first  introduced  in 
[3].  The  implementation  of  this  idea  for  a  chain  of 
spins  interacting  through  the  Ising  interaction  was 
given  in  [4].  Ising-type  spin  chains  have  been  used 
for  quantum  computation  in  a  statistical  ensemble  of 
quantum  computers  [5,6].  In  [7],  we  considered  a  sim¬ 
ple  model — a  chain  of  identical  spins  in  a  non-uni  form 
magnetic  field. 

The  present  paper  is  based  on  [4,7].  We  con¬ 
sider  a  chain  of  spins  (e.g.,  nuclear  spins)  in  a  non- 
uniform  magnetic  field.  Similar  to  [7],  the  angle,  (9, 
between  the  direction  of  the  chain  and  the  direction 
of  the  permanent  magnetic  field  (z-dircction)  satis¬ 
fies  the  condition:  cos  (9  =  l/\/3.  Then,  the  dipole- 
dipole  interaction  between  spins  is  suppressed,  and  the 
Ising  interaction  becomes  dominant.  We  suppose  that 
our  chain  consists  of  a  periodic  sequence  of  triplets, 
ABC  ABC  ABC  —  The  triplet,  ABC ,  can  be  differ¬ 
ent  nuclear  spins  or  identical  spins  in  slightly  different 
environments.  For  definiteness,  we  shall  keep  in  mind 


326 


G.P  Berman  et  al.  /  Computer  Physics  Communications  146  (2002)  324-330 


this  second  case  [7]  with  the  following  typical  para¬ 
meters, 

co\  /2n  ~  400  MHz,  co*  =  co\  +  (k  —  1)Acd, 

Au)/2tt  %  20  kHz,  (5) 

Jac/2jtk\00  Hz,  Jbc  =  2Jac ,  Jab  =  2>Jac > 

where  co\  is  the  Larmor  (angular)  frequency  of  the 
right  end  spin,  a;*  is  the  Larmor  frequency  of  the  fcth 
spin,  Jik  is  the  constant  of  Ising  interaction  between 
the  / th  and  kth  spins. 

In  the  presence  of  a  circularly  polarized  transverse 
magnetic  field,  the  Hamiltonian  of  the  system  can  be 
written  as  [8], 

n  =  -J2MkIk  ~2Y1  y*  *+i 7*  7*+t 

k  k 

-  y  J2 1 4~  exP[-'M  +  <*>)] 

k 

+  /*+exp[i(ajf  +  <p)]},  (6) 

where  lz  is  the  nuclear  spin  operator,  co  and  tp  are 
the  frequency  and  the  phase  of  the  transverse  mag¬ 
netic  field;  £2  is  the  Rabi  frequency  (the  magnitude  of 
the  transverse  magnetic  field  in  frequency  units),  and 
h  =  1. 

5.  Implementation  of  the  quantum  full  adder 

First  we  shall  define  our  problem.  Suppose  that 
we  have  a  number  | b)  =  \b^L~^  . . .  M0)),  b =0,1. 
In  decimal  notation,  b  =  b^L~l)2L~]  +  •  •  •  +  b(0) 2°. 
(Below,  we  omit  parentheses  in  the  superscripts.) 
Note,  that  in  general,  we  have  a  superposition  of  many 
numbers,  \b)i ,  and  any  gate  must  act  on  all  of  these 
numbers  simultaneously.  We  are  going  to  add  to  a 
number  | b)  (all  numbers,  bi ,  in  a  superposition),  a 
definite  number,  a  =  aL~x  . .  .a0,  where  am  =  0,  1 .  To 
achieve  this  goal  we  shall  use  2 L  +  1  qubits.  We  load 
the  number  b  in  a  chain  of  qubits  in  the  following  way, 

\02L+ib^['02L-\b^2 . .  .05fc]0302^)-  (7) 

This  means  that  we  place  two  additional  qubits  in  the 
states  |0)  in  front  of  the  qubit  ,  and  one  additional 
qubit  in  the  state  |0)  in  front  of  all  other  qubits  bm , 
/n/0. 

The  gate  C 123  is  not  convenient  for  our  spin  chain 
quantum  computer  in  which  each  spin  interacts  only 


with  its  neighbors.  So,  we  replace  it  with  the  following 
sequence  of  gates, 

c  123  =  C23C32C23C132C23C32C23.  (8) 

Let  us  check  this  equation,  for  example,  for  the  state, 
IO3I21 1).  From  the  left  side  of  Eq.  (8)  we  have, 

123IO3 12 1 1  >  =  1 13 12 1 1 ),  (9) 

as  the  first  and  the  second  qubits  are  in  the  states,  1 1), 
the  third  qubit  changes  its  state.  Next,  we  follow  the 
operations  on  the  right-hand  side  of  (8)  to  show  that 
the  same  result  is  obtained.  From  the  right  side  of 
Eq.  (8)  we  have, 

C23|03l2ll>  =  |l3l2ll>,  C321 1 3 1 2 1 1  >  =  1 1 3O2 1 1 )  * 

<^23 1 1 3O2  1 1 )  =  1 1 3O2  1 1  >,  C*  I32l  I3O2 1 1  >  =  1 13 12  1 1  >• 

C23II3I2I1)  =  IO3  12  1 1  >•  C32IO3I2I1)  =  IO3I2I1), 

C23IO3I2I1)  =  1 13  12  1 1 ), 

(10) 

which  coincides  with  the  right  side  of  Eq.  (9).  Now, 
instead  of  (2b),  we  have  the  following  expression  for 
the  full  adder,  F, 

F321  (0)  =  C12C23C32C23C132C23C32C23, 

^32l(l)  =  C'i2C23C32C23Ci32C23C32C23C,2C23.  ^ 

Next,  we  explain  how  to  add  the  numbers  b  and 
a.  If  a0  =  0,  the  quantum  computer  applies  the  gate, 
F321  (0).  If  a0  =  1,  the  quantum  computer  applies  the 
gate,  F321O).  According  to  (2),  the  result  is, 

^32i|0302*?)  =  |c]jX).  (12) 

where  we  have  replaced  s'  by  s°  =  a0  0  b°y  and  c'  by 
c]  which  is  the  carry-over  for  the  next  addition  of  b] 
and  a 1 . 

Next,  consider  the  five  right  most  qubits, 

IO5  b\c\s^b\).  (13) 

To  add  bx  and  ax  and  the  carry-over,  c1,  we  should 
first  make  a  swap,  5,  of  the  values  of  the  fourth  and 
the  third  qubits, 

Su\b\c\)  =  \c\b\).  (14) 

The  swap  gate,  S,*,  can  be  represented  in  terms  of  C,* 
gates, 

Sik  =  CkiCikCki. 


(15) 
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Let  us  check,  for  example,  the  action  of  the  521  gate 
on  the  state,  IO2 1 1 )-  Using  (15),  we  have, 

Si\  =  C12C21C12, 

Ci2|02li>  =  |l2li),  U2i  1 1 2 1 1  >  =  1 1 20 1  > ,  (16) 

C12II2OO  =  |l20i>. 

Thus,  the  5-gate  transforms  the  state,  |01),  into  the 
state,  1 10). 

After  the  action  of  the  543  gate,  the  state  (13) 
changes  into, 

|0  5c'4b\sffi).  (17) 

The  state,  IO5 c\b\),  has  the  form  (2a),  and  it  is  ready 
for  application  of  the  full  adder,  F543, 

Fsu\05c\bl)  =  \c25s\b\).  (18) 

Certainly,  we  use  ^543(0)  if  a1  =  0  and  F543O)  if 
ax  =  1.  Thus,  we  obtain  the  sum,  s1  =  a1  ©  b]  0  c1, 
and  a  carry-over,  c2,  for  the  next  addition.  It  is  clear 
that  by  repeating  the  application  of  gates  F  and  5, 
we  shall  obtain  the  desired  answer.  A  complete  full 
adder,  F,  can  be  represented  as  a  combination  of  the 
elementary  full  adders,  Fy*,  and  the  swap-gates,  5y, 

F  =  F2L+\%2L,2L-\(aL~X)S2L2L-\  •  •  • 

F765  (a 2 )  $65  F543  (a 1 )  543  F32 1  (a 0 ) .  (19) 

After  the  action  of  the  full  adder,  the  superposition 
of  the  states  (7)  transforms  into  the  superposition  of 
states, 

I  cLsL-'bL-'sL-2bL-2...s°b0).  (20) 

Thus,  the  qubits  in  even  positions  and  the  left  end 
qubit  carry  the  results  of  addition.  The  problem  is: 
How  to  implement  all  of  these  gates  using  ;r -pulses? 
We  suppose  that  Q  <£  7/*,  say  £2/2: r  ^  10  Hz. 
This  allows  us  to  excite  each  spin  individually.  To 
implement  F321,  we  use  Eqs.  (11).  As  an  example, 
we  show  how  to  implement  F32i(0).  According  to 
Eqs.  (11),  first  we  implement  C23,  then  C32,  and  so 
on.  Suppose  we  have  an  integer  number  of  triples, 
ABC ,  in  our  spin  chain.  Then,  we  have  4  possible 
frequencies  in  the  position  (3 n  —  2)  which  correspond 
to  a  spin  C: 

&>3^_2  =  <^3rt-2  4-  JbC  4-  j AC, 

"3«-2="3n— 2  +  JBC  ~  JAC,  (2j) 

tt>3®_2  =  <*>3n-2  -  JBC  +  J AC . 

=&>3n-2  -  JbC  ~  JaC ■ 


Here  col£  corresponds  to  the  k\h  spin  whose  left 
neighbor  is  in  state  | i)  and  whose  the  right  neighbor  is 
in  state  \j ).  Similar  expressions  can  be  found  for  spins 
in  positions  (3 n  —  1)  and  (3 n),  the  B  and  A  spins, 
respectively.  The  end  spins  have  only  two  frequencies. 
For  the  left  end  spin,  A, 

^2L+i  =  ^2 l + 1  +  Jab,  ^2l+i  =Cl)2c+i  ~  Jab, 

(22) 

and  for  the  right  end  spin,  C, 

co®  =  co\  +  J  bc,  (o\  =  co\  —  J  bc-  (23) 

Now,  to  implement  C23,  we  apply  a  7r -pulse  with 
frequency  co®1  and  then  a  7r -pulse  with  the  frequency 
cu]1.  One  of  these  two  pulses  definitely  changes  the 
state  of  the  third  spin  if  the  second  spin  is  in  its  excited 
state,  |1).  To  implement  the  gate  C123,  we  need  a 
single  71  -pulse  with  frequency  cu]1.  To  implement  a 
Not-gate,  C2,  which  appears  in  F132U),  we  have  to 
apply  four  ;r -pulses  with  all  possible  frequencies,  co ^  , 
i  =  0,  1,  and  j  =0,1.  The  total  number  of  pulses 
required  to  implement  a  F/;*(0)-gate  (if  /  ^  2L  +  1, 
k  7^  1)  is  15:  two  pulses  for  each  Control-Not  gate 
and  one  pulse  for  the  Control-Control-Not  gate.  For 
the  Fjjid  1)  gate,  the  total  number  of  pulses  is  21  (16 
pulses  for  the  left  triple  of  the  chain).  The  swap  gate, 
Sij ,  requires  6  pulses. 

Thus,  to  add  L-qubit  numbers,  our  proposal  re¬ 
quires  (2 L  +  1)  qubits  and  less  than  27 L  n -pulses. 

6.  Simulations  of  a  quantum  full  adder 

For  numerical  simulations  of  a  quantum  full  adder 
we  used  the  following  assumption:  the  frequency 
difference  between  two  neighboring  spins  is  much 
greater  than  the  Rabi  frequency.  As  a  result,  the 
selective  excitation  of  a  chosen  spin  has  a  small  effect 
on  all  other  spins.  At  the  same  time,  we  take  into 
account  the  action  of  a  tt  -pulse  on  non-resonant  states. 
For  example,  suppose  that  the  frequency  of  a  n -pulse 
is:  co  =  co\x .  This  pulse  is  resonant  with  spin  “3”  only 
if  the  neighboring  spins,  “1”  and  “2”,  are  in  their 
excited  states.  The  states  in  which  the  spin  “3”  has 
frequencies,  cu]°,  co® 1 ,  and  CO30,  are  non-resonant  for 
the  pulse  with  co  =  co\ 1 .  We  take  into  consideration  the 
transformations  of  all  these  non-resonant  states. 
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Note,  that  any  pulse  in  our  simulations  acts  on 
all  basic  states  in  the  quantum  superposition.  The 
resonant  state  transforms  into  a  state  with  the  opposite 
direction  of  the  resonant  spin.  Every  non-resonant 
state  transforms  into  the  superposition  of  two  states: 
the  initial  one  and  an  error  state  generated  by  the  non¬ 
resonant  transition.  If  the  probability  of  an  error  state 
is  less  than  a  chosen  small  number,  e ,  our  computer 
program  automatically  removes  it.  One  can  argue  that 
every  removed  state  can  generate  a  number  of  new 
states.  However,  suppose  that  we  have  removed  a 
state  with  a  small  probability,  P.  Due  to  the  main 
property  of  the  unitary  transformations,  the  sum  of 
the  probabilities  of  all  states  generated  by  the  removed 
state  (including  itself)  is  P.  Thus,  the  total  probability 
of  all  neglected  states  cannot  increase  in  spite  of 
generation  of  new  states! 

We  have  simulated  the  addition  of  the  100-digit 
numbers  in  a  201  spin  chain.  To  minimize  non¬ 
resonant  errors,  we  used  the  27r£;-method  [7,8].  The 
basic  idea  of  the  27r/:-method  is  the  following.  One 
chooses  the  Rabi  frequency,  f2,  of  a  resonant  tt -pulse 
in  such  a  way  that  it  becomes  approximately  a  2i re¬ 
pulse  for  all  non-resonant  transitions  (where  k  is  an  in¬ 
teger  which  generally  is  different  for  different  states). 
A  27r&-pulse  does  not  generate  unwanted  error  states. 
A  transformation  of  a  basic  state  under  the  action  of  a 
7r-pulse  was  described  in  our  previous  paper  [7].  The 
following  values  of  dimensionless  parameters  were 
chosen, 

Jac  =  1,  Jbc  =  2,  JAB  =3.  (24) 

We  found  that  all  non-resonant  transitions  approxi¬ 
mately  satisfied  the  27r/:-conditions  for  the  value  of 
the  Rabi  frequency,  Q  =  Qo  =  0. 10005. 

Next,  we  present  the  results  of  our  numerical 
simulations.  As  an  example,  we  add  the  “classical 
number”, 

199098.. .0°  (25) 

(299  in  a  decimal  notation)  to  the  “quantum  number”, 

O99 ...  0 1 1 0  (26) 

(1  in  a  decimal  notation).  The  sum  of  these  two 
numbers, 

199098 . .  .01 1°,  (27) 

corresponds  to  the  quantum  state  of  the  chain  of  201 
spins, 


IO201 1200O199 . . .  O3 12 1 1  >-  (28) 

For  the  value  of  =  £2q  which  approximately 
satisfies  27r/:-condition,  we  have  for  the  probability  of 
the  expected  state  (28):  Po  =  0.99889.  The  number, 
N ,  of  error  states  with  probability  P  >  10“ 12  is 
only  304.  Small  deviations  from  the  27r/:-condition 
significantly  influence  the  result.  As  an  example,  for 
f2  =  0.10021,  the  value  of  Po  decreases  to  0.98300, 
and  the  number  of  error  states,  A,  grows  to  46530. 
Fig.  1  shows  the  dependence  of  Po  and  N  on  the 
Rabi  frequency,  in  the  vicinity  of  f2o-  One  can 
see  that  the  probability,  Po ,  smoothly  decreases  by 
approximately  1%  when  the  Rabi  frequency,  f2,  shifts 
from  S2q  by  approximately  0.1%.  Unlike  the  value 
of  P0,  the  number  of  error  states  with  probability 
>  10~12  does  not  change  in  the  close  vicinity  of  &q. 
But  it  sharply  increases  when  the  deviation  of  J2  from 
f2o  approaches  0.15%.  Fig.  2  shows  the  probabilities 


Fig.  1.  The  probability  of  the  expected  state,  Pq  ;  and  the  number  of 
error  states,  N,  as  a  function  of  the  Rabi  frequency,  C2 . 
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Fig.  2.  “Line  spectrum”  of  the  probabilities,  P,  of  error  states 
(f2  =  0.10021);  (a)  the  region:  P  ~  10-6,  (b)  the  region: 
P~  10"10. 
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of  error  states  (the  states  are  shown  in  the  order  of  their 
generation).  One  can  see  a  specific  “line  spectrum” 
of  the  probabilities:  the  error  states  are  “attracted”  to 
a  few  discrete  values  of  the  probability.  Very  similar 
effects  were  obtained  for  other  examples  of  addition 
including  the  addition  of  a  “classical”  number  to  a 
superposition  of  two  quantum  numbers. 


7.  Phase  control 

The  quantum  full  adder  implementation  considered 
above  provides  a  transformation  of  an  arbitrary  super¬ 
position  of  “quantum  numbers”,  qk ,  into  a  superposi¬ 
tion  of  numbers  (qk  +  a),  where  a  is  any  given  “clas¬ 
sical  number”.  However,  this  proposed  scheme  also 
generates  complicated  phase  differences  between  the 
states  of  the  quantum  superposition.  This  effect  can 
be  inappropriate,  especially  for  the  Shor  algorithm.  In 
this  section,  we  describe  how  to  extend  our  simple 
model  to  incorporate  phase  restoration  after  the  action 
of  every  tt  -pulse. 

Consider  the  chain  of  paramagnetic  ions  ABCABC 
containing  nuclear  spins  in  a  high  nonuniform  mag¬ 
netic  field.  Every  electron  spin  interacts  with  the  nu¬ 
clear  spin  of  its  ion  via  the  hyperfine  interaction.  Elec¬ 
tron  spins  interact  with  each  other  through  the  Ising  in¬ 
teraction.  All  electron  spins  are  in  their  ground  states. 
The  nuclear  spins  also  interact  with  each  other  through 
the  Ising  interaction.  This  interaction  is  responsible  for 
quantum  logic  gates.  The  dipole-dipole  interaction  is 
suppressed  due  to  the  use  of  the  magic  angle  between 
the  chain  and  the  direction  of  the  external  magnetic 
field.  The  key  point  of  our  model  is  the  following:  we 
consider  the  interaction  (for  simplicity,  the  Ising  inter¬ 
action)  between  an  electron  spin  and  its  two  neighbor¬ 
ing  nuclear  spins.  This  interaction  can  originate,  for 
example,  if  the  electron  density  at  the  neighboring  nu¬ 
clei  is  not  zero.  Thus,  the  electron  spin  frequency  for 
a  particular  ion  can  take  eight  values.  As  an  example, 
for  the  B-ion  one  has, 

CO000  =a>e  +  JeAeB  +  JBC  +  AB/2  +  JAB  +  JBC, 
w100  =coe  +  JAB  +  JBC  -  AB/2  +  JAB  +  JeBC , 
a,1 10  =  a*  +  JAB  +  JBC  -Ab/2-  Jab  -  JBC, 

(29) 


and  so  on.  Here  a/7*  corresponds  to  the  electron 
frequency  for  the  case  in  which  the  nuclear  spin  of 
the  same  ion  is  in  the  state  |i),  the  nuclear  spin  of 
the  neighboring  A-ion  is  in  the  state  | k),  and  the 
nuclear  spin  of  the  neighboring  C-ion  is  in  the  state 
|  j)\  coe  is  the  Zeeman  frequency,  J^B  and  JBec  are 
the  constants  of  the  electron-electron  interaction,  AB 
is  the  hyperfine  constant  (in  the  frequency  units),  and 
J^B  and  JBnc  are  the  constants  of  the  electron-nuclear 
interaction  for  neighboring  ions. 

Thus,  the  ESR  frequency  depends  on  the  position 
of  an  electron  spin  in  the  chain  (because  the  Zeeman 
frequency  is  nonuniform)  and  the  states  of  three 
nuclear  spins — the  nuclear  spin  of  its  own  ion  (via 
the  hyperfine  interaction)  and  the  nuclear  spins  of  two 
neighboring  ions.  One  can  tune  the  frequency  of  an 
electromagnetic  pulse  in  such  a  way  that  it  is  resonant 
with  the  electron  spin  only  if  it  is  in  a  definite  position 
in  the  chain  and  the  three  nuclear  spins  mentioned 
above  are  in  definite  states. 

The  strategy  for  phase  correction  for  the  quantum 
full  adder  is  the  following.  The  full  adder  is  imple¬ 
mented  by  a  sequence  of  the  nuclear  tt  -pulses.  A  nu¬ 
clear  tt -pulse  causes  a  phase  shift  of  tt/2  for  the  reso¬ 
nant  states.  There  are  six  possible  phase  distortions  for 
non-resonant  states  (3  possible  states  for  neighboring 
nuclear  spins  times  two  possible  states  for  the  selected 
nuclear  spin).  Correspondingly,  one  has  six  possible 
frequencies  for  the  electron  spin  of  the  selected  ion  in 
the  non-resonant  state.  Because  of  the  large  value  of 
the  electron  gyromagnetic  ratio,  we  assume  that  each 
of  the  corresponding  electron  transitions  can  be  driven 
without  noticeable  non-resonant  effects  on  the  elec¬ 
tron  transitions  with  close  frequencies.  After  the  nu¬ 
clear  7i -pulse,  one  applies  12  electron  tt -pulses:  two 
7i  -pulses  for  every  possible  frequency  of  the  selected 
electron  spin  in  the  non-resonant  state.  The  phase  of 
the  first  7i  -pulse  is  zero,  the  phase  of  the  second  one 
is  0.  The  total  phase  shift  of  the  wave  function  of  the 
ion  is  (tt  +  0).  By  choosing  an  appropriate  value  for 
0,  one  can  change  the  phase  shift  for  a  specific  non¬ 
resonant  state  to  the  value  of  7r/2,  which  corresponds 
to  the  resonant  state.  After  the  action  of  12  electron 
7i -pulses,  all  states  in  the  quantum  superposition  will 
have  the  same  phase,  and  all  electron  spins  will  be  in 
their  ground  state.  Thus,  the  phase  distortion  generated 
by  a  nuclear  tt  -pulse  can  be  corrected  with  1 2  electron 
tt  -pulses. 
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8.  Conclusions 

We  have  demonstrated  a  strategy  for  implementa¬ 
tion  of  a  full  adder  in  a  nuclear  spin  quantum  com¬ 
puter.  To  add  an  arbitrary  superposition  of  L-qubit 
“quantum  numbers”  and  a  fixed  L-bit  “classical  num¬ 
ber”,  our  scheme  requires  2L  -f  1  qubits  and  less  than 
27 L  resonant  tt  -pulses.  We  have  simulated  the  ac¬ 
tion  of  this  full  adder  for  a  spin  chain  containing  201 
qubits.  Using  the  27r/:-method  allowed  us  to  minimize 
the  generation  of  unwanted  non-resonant  error  states. 
These  error  states  are  shown  to  have  a  few  preferred 
discrete  values  of  probability  (a  “line  spectrum”).  Our 
simple  spin  model  is  generalized  to  correct  undesired 
phase  distortions. 
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Abstract 

We  survey  recent  work  on  designing  and  evaluating  quantum  computing  implementations  based  on  nuclear  or  bound-electron 
spins  in  semiconductor  heterostructures  at  low  temperatures  and  in  high  magnetic  fields.  General  overview  is  followed  by  a 
summary  of  results  of  our  theoretical  calculations  of  decoherence  time  scales  and  spin-spin  interactions.  The  latter  were  carried 
out  for  systems  for  which  the  two-dimensional  electron  gas  provides  the  dominant  carrier  for  spin  dynamics  via  exchange  of 
spin-excitons  in  the  integer  quantum  Hall  regime.  ©  2002  Elsevier  Science  B.V.  All  rights  reserved. 

PACS:  73.20.Dx;  71.70.Ej;  03.67.Lx;  76.60.-k 
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1.  Introduction 

The  field  of  quantum  computing  has  seen  explo¬ 
sive  growth  of  experimental  and  theoretical  interest. 
The  promise  of  quantum  computing  [1-5]  has  been  in 
exponential  speedup  of  certain  calculations  via  quan¬ 
tum  parallelism.  In  Fig.  1,  the  top  flow  chart  shows 
the  “classical”  computation  which  starts  from  binary 
input  states  and  results  in  binary  output  states.  The  ac¬ 
tual  dynamics  is  not  really  that  of  Newtonian  classi¬ 
cal  mechanics.  Rather  the  computation  involves  many- 
body  irreversible  “gate”  device  components,  made  of 
semiconductor  materials  in  modem  computers,  which 
evolve  irreversibly,  “thermodynamically”  according  to 
the  laws  of  statistical  mechanics.  As  the  size  of  the 
modem  computer  components  approaches  atomic,  the 
many-body  quantum  behavior  will  have  to  be  ac¬ 
counted  for  in  any  case  [6]. 
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The  idea  of  quantum  computing,  however,  is  not 
just  to  account  for,  but  to  actually  utilize  the  quantum- 
mechanical  dynamical  behavior.  This  is  not  an  easy 
task.  Quantum  mechanics  allows  for  parallelism  in 
evolution:  one  can  “process”  a  linear  superposition 
of  several  input  states  at  once,  as  illustrated  in  the 
lower  flow  chart  in  Fig.  1.  The  price  paid  is  that 
coherent  processing  of  information,  according  to  the 
law  of  quantum  mechanics,  must  be  accomplished 
in  systems  much  larger  than  atomic-size  (or  more 
importantly,  with  many  degrees  of  freedom).  There 
arc  numerous  conceptual  and  experimental  obstacles 
to  accomplishing  this  task,  that  have  generated  a  lot 
of  interest,  excitement,  and  new  results  in  computer 
science,  physics,  and  engineering. 

The  functioning  of  a  quantum  computer  involves 
initialization  of  the  input  state,  then  the  actual  dy¬ 
namical  evolution  corresponding  to  computation,  and 
finally  reading  off  the  result.  Various  specific  re¬ 
quirements  for  implementation  have  been  identified 
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Fig.  1 .  Comparison  of  the  classical  and  quantum  approaches  to  com¬ 
puting.  The  upper  flow  chart  schematically  represents  implemen¬ 
tation  of  a  traditional  irreversible  “classical”  computation  process, 
where  transformation  of  the  input  set  of  bits  into  the  result  is  ac¬ 
complished  by  a  succession  of  irreversible  gates.  Owing  to  their  irre¬ 
versibility,  the  gates  can  be  connected  in  space  rather  than  switched 
on  and  off  at  different  times.  The  lower  flow  chart  shows  quantum 
processing  of  information,  where  the  input  and  the  final  result  are 
both  in  superposition  states,  yielding  quantum  parallelism.  The  dy¬ 
namics  is  reversible:  there  is  a  one-to-one  correspondence  between 
the  initial  and  final  states.  Therefore,  number  of  the  input  and  out¬ 
put  quantum  bit  (qubits)  is  the  same  even  though  some  of  the  output 
qubits  (set  in  a  smaller  font)  might  not  be  used  in  the  final  extrac¬ 
tion  of  the  classical  result  by  measurement.  The  quantum  gates  are 
applied  in  succession  by  being  switched  on  and  of  at  different  times 
during  the  computation.  The  question  mark  signifies  the  difficulty  of 
finding  quantum  algorithms  that  retain  the  power  of  quantum  paral¬ 
lelism  after  measurement  needed  to  read  off  the  final  result  as  clas¬ 
sical  information. 


[2-5];  here  we  provide  only  a  limited  introductory 
overview. 

Let  us  begin  by  considering  the  reading  off  of  the 
final  result.  The  reason  for  the  question  mark  in  the 
lower  chart  in  Fig.  1  is  that  quantum  measurement  of 
the  final  superposition  state  can  erase  the  gain  of  the 
parallel  dynamics,  by  collapsing  the  wave  function. 
Therefore,  a  key  issue  in  quantum  computing  has 
been  to  find  those  algorithms  for  which  the  readout 
of  the  final  state,  by  way  of  projecting  out  a  certain 
average  property,  still  retains  the  power  of  the  quantum 
parallelism.  To  date,  only  few  such  examples  are 
known  [1,3, 4, 7],  the  most  celebrated  being  the  Shor 
algorithm  [  1  ]  for  factoring  of  integers,  the  invention  of 
which  boosted  quantum  computing  from  an  obscure 
theoretical  field  to  a  mainstream  research  topic. 

The  preparation  of  the  initial  state  does  not  seem 
to  present  a  problem  for  most  quantum  computing 
realizations  [2-5],  except  perhaps  the  ensemble  liq¬ 


uid  state  NMR  approach  [8,9]  which  relies  on  the 
initial  thermal  distribution  to  produce  deviation  of 
the  density  matrix  from  the  equal-probability  mixture 
state.  In  most  other  approaches,  the  initial  state  can 
be  produced  by  first  fully  polarizing  the  quantum  bits 
(qubits),  i.e.  putting  them  in  one  of  the  two  quantum 
levels.  Note  that  we  consider  two-state  qubits  here,  re¬ 
alized,  for  instance,  by  spins  1/2  of  nuclei  or  gate- 
or  impurity-bound  electrons,  in  applied  magnetic  field. 
The  fully  polarized  state  is  then  subject  to  gate  opera¬ 
tions  to  form  the  desired  input  state.  Part  of  a  quantum¬ 
computing  algorithm  should  be  the  prescription  on 
how  to  choose  the  initial  state  to  represent  the  clas¬ 
sical  information  of  the  input,  like  the  input  integer  in 
the  factoring.  In  most  cases,  this  prescription  is  easily 
accomplished  by  single-qubit  and  two-qubit  gates. 

The  actual  dynamical  evolution  (the  process  of 
computation)  in  quantum  computing  is  fully  reversible 
and  nondissipative,  unlike  classical  computing.  Much 
progress  has  been  made  in  resolving  both  the  con¬ 
ceptual  and  computer-engineering  “design”  issues  for 
quantum  computation.  Specifically,  the  computation 
can  be  carried  out  [2-5,10-13]  by  a  universal  set  of 
gates:  single-qubit  rotations  and  nearly  any  two-qubit 
gate.  The  gates  are  not  connected  in  space  like  in  clas¬ 
sical  computers  but  are  activated  in  succession  in  time, 
to  control  single-spin  dynamics  and  also  switch  on  and 
off  two-spin  interactions  (we  use  “spin”  and  “qubit” 
interchangeably). 

Many  interesting  matters  have  been  resolved,  which 
are  not  reviewed  here.  These  include  the  understand¬ 
ing  of  how  the  finiteness  of  the  state  space  (i.e.  two 
states  for  spin  1/2)  replaces  the  “classical”  digitaliza¬ 
tion  in  quantum  computing.  Also,  the  “classical”  copy¬ 
ing  (fan-out)  function  is  not  possible  in  quantum  me¬ 
chanics.  It  is  replaced  by  entanglement  with  ancillary 
qubits  to  accomplish  redundancy  needed  for  error  cor¬ 
rection  [14-20].  Sources  of  errors  due  to  interactions 
with  environment  in  quantum  mechanics  involve  not 
only  the  usual  relaxation  (thermalization)  but  also  loss 
of  coherence  [21-28].  This  quantum  decoherence  (de¬ 
phasing)  can  be  faster  than  relaxation  because  it  does 
not  require  energy  exchange. 

A  conceptually  important  issue  has  been  the  scala¬ 
bility  of  quantum  computing:  can  one  process  macro- 
scopically  large  amounts  of  information  by  utilizing 
quantum  error  correction  based  on  redundancy  via  en¬ 
tanglement  with  ancillary  qubits?  The  affirmative  an- 
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swer  to  this  question  has  been  one  of  the  triumphs  of 
the  theory  [14-20].  It  provided  a  new  paradigm  for 
emergence  of  controlled/organized  macroscopic  be¬ 
havior  from  microscopic  dynamics,  on  par  with  the 
conceptual  possibility  of  living  organisms,  which  we 
observe  but  cannot  yet  “manufacture”,  and  million- 
gate  classical  computers  which  are  man-made. 

With  all  these  theoretical  advances  at  hand,  the  next 
step  is  to  ask  whether  a  man-made  quantum  computer 
can  be  realized?  There  have  been  several  experimental 
directions  of  exploration,  most  presently  are  still  at  the 
level  of  one  or  two  qubits,  or,  for  ensemble  liquid- 
state  NMR,  which  emulates  quantum  dynamics  by 
evolution  of  the  density  matrix  of  a  large  collection 
of  molecules,  5-7  qubits. 

In  this  introductory  survey,  we  summarize  results 
of  our  work  on  two-spin  interactions  and  spin  decoher¬ 
ence  in  semiconductor  heterostructures.  In  Section  2, 
we  consider  the  spin-based  quantum  computing  pro¬ 
posals  in  such  systems.  Time  scales  of  relaxation  and 
decoherence  are  addressed  in  Section  3.  Finally,  Sec¬ 
tion  4  reports  results  for  models  with  nuclear  spins  as 
qubits. 

2.  Spin-based  quantum  computing  in 
semiconductor  heterostructures 

The  general  layout  of  a  solid-state  quantum  com¬ 
puter  is  shown  in  Fig.  2.  Qubits  are  positioned  with 
precision  of  few  nanometers  in  a  heterostructure.  One 
must  propose  how  to  effect  and  control  single-qubit 
interactions,  two-qubit  interactions,  and  explore  how 
the  controlled  dynamics  owing  to  these  interactions 
compares  to  decoherence  and  relaxation.  The  proposal 
must  include  ideas  for  implementation  of  initializa¬ 
tion,  readout,  and  gate  functions. 

The  first  proposal  including  all  these  components 
was  for  qubits  realized  in  an  array  of  quantum  dots 
[29]  coupled  by  electron  tunneling.  The  first  spin- 
based  proposal  [30]  utilized  nuclear  spins  coupled  by 
the  two-dimensional  electron  gas,  the  latter  in  the  dis¬ 
sipationless  integer  quantum  Hall  state  [31]  that  re¬ 
quires  low  temperatures  and  high  magnetic  fields.  An 
important  advancement  was  the  work  of  Kane  [32] 
where  gate  control  of  nuclear-spins  of  donor  impuri¬ 
ties,  separated  less  than  10  nm  and  coupled  via  the 
outer  impurity  electrons  which  are  bound  at  low  tem- 


Fig.  2.  Schematic  illustration  of  a  semiconductor  heterostructure 
quantum  information  processor.  The  qubits,  represented  by  the 
arrows  overlaying  heavy  dots,  are  spins  1/2  of  nuclei  or  localized 
electrons.  Individual  control  of  the  temporal  evolution  of  the  spins 
can  be  achieved  with  the  use  of  external  electromagnetic  radiation, 
i.e.  NMR  or  ESR  pulses.  The  spins  are  also  coupled  with  each 
other  via  interaction  mediated  by  the  two-dimensional  electron 
gas  in  the  heterostructure,  or  by  other  means.  The  external  and 
internal  interactions  can  be  controlled  by  gates  formed  on  top  of 
the  heterostructure.  The  external  environment,  that  includes  crystal 
lattice,  electron  gas,  defects,  impurity  potentials,  causes  relaxation 
and  decoherence  of  the  qubits. 

peratures,  was  proposed.  Most  of  these  ideas  also 
apply  to  electron-spin  qubits,  bound  at  impurities, 
in  quantum  dots,  or  directly  by  gates.  Several  elab¬ 
orate  solid-state  heterostructure  quantum  computing 
schemes  have  been  proposed  in  the  literature  recently 
[28,33-41].  There  are  also  other  promising  proposals 
involving  surface  geometries:  superconducting  elec¬ 
tronics  [42—46]  and  electrons  on  the  surface  of  liquid 
helium  [47]. 

There  have  been  several  planned  and  ongoing 
experimental  efforts  [32-36,43-45,48-54]  ultimately 
aimed  at  solid-state  quantum  computing  and  other 
quantum  information  processing  realizations.  The  fi¬ 
nal  geometry  is  expected  to  be  most  sensitive  to  the 
implementation  of  readout,  because  it  involves  quan¬ 
tum  measurement,  i.e.  supposedly  interaction  with 
or  transfer  of  information  to  a  macroscopic  device. 
Therefore,  much  of  the  experimental  effort  presently 
has  been  focused  on  single-qubit  (single-spin)  mea¬ 
surement  approaches. 

The  theoretical  efforts  can  be  divided  into  two  ma¬ 
jors  tasks.  The  process  of  single-spin  measurement 
must  be  understood  for  the  readout  stage  of  quantum 
computing.  Several  conceptual  and  calculational  ad¬ 
vancements  have  been  made  in  understanding  quan¬ 
tum  measurement  [26,32-36,46,50,55,56]  as  it  applies 
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to  atomic-size  qubit  systems  interacting  with  environ¬ 
ment  and  typically  “measured”  directly  by  the  effect 
of  the  spin-qubit  state  on  transport,  or  first  transferring 
the  spin  state  to  a  charge  state  that  is  easier  to  measure, 
e.g.,  in  single-electron  transistors  and  similar  devices. 

In  this  survey,  we  outline  results  of  the  second 
evaluation  task:  that  of  understanding  the  processes 
and  times  scales  involved  in  the  dynamics  of  the  actual 
computation.  As  summarized  in  Fig.  3,  this  main 
stage  of  the  quantum  computation  process  involves 
control  of  spins  and  their  interactions.  It  also  involves 
processes  that  we  do  not  control  and  are  trying  to 
minimize:  relaxation  and  decoherence. 

Control  of  individual  qubits  is  usually  accom¬ 
plished  externally.  For  nuclear  spins,  NMR  radio¬ 
frequency  radiation  can  be  used,  see  Fig.  2.  For  elec¬ 
tron  spins,  the  ESR  microwave  frequencies  are  suit¬ 
able.  Such  radiation  cannot  be  focused  on  the  scale 
of  10-100  nm.  Instead,  selectivity  must  be  accom¬ 
plished  by  independent  means.  Several  proposals  ex¬ 
ist,  the  most  promising  being  control  by  gates.  The  ap¬ 
plied  gate  voltage  modifies  the  electronic  wave  func¬ 
tion  changing  interactions  and  therefore  resonant  fre¬ 
quencies.  We  will  denote  the  time  scale  of  the  external 
single-qubit  control  by  Tex{.  This  can  be  the  Rabi  time 
of  a  spin  flip. 

The  qubit-qubit  interactions  are  typically  assumed 
to  be  mediated  by  electrons  that  “visit”  both  qubit 
environments.  For  instance,  in  liquid-state  ensemble 
NMR  [8,9]  with  complex  molecules,  or  in  the  original 
model  [32]  of  phosphorous  impurity  donors  in  silicon, 
the  wave  functions  of  the  valence,  outer  electrons  of 
nearby  qubits  overlap.  Specifically,  in  the  P  donor 
case,  the  single  outer  electron  of  the  donor  atom 
remains  bound  at  low  temperatures  but  has  orbital 
radius  of  order  2  nm  owing  to  the  large  dielectric 
constant  of  the  silicon  host.  Therefore,  it  is  hoped 
that  these  electrons,  in  nearby  donors  positioned  as  in 
Fig.  2,  will  mediate  nuclear-spin  qubit  interactions. 

Our  approach  [27,28]  allows  for  larger  qubit  sep¬ 
aration,  up  to  order  100  nm,  by  relying  on  the  two- 
dimensional  electron  gas  in  the  heterostructure  to  me¬ 
diate  qubit-qubit  interactions.  This  two-dimensional 
electron  gas  is  usually  obtained  by  spontaneous  or 
gate-induced  transfer  of  electrons  from  impurities 
to  the  two-dimensional  interface  layer  in  which  the 
qubits  are  positioned.  The  source  impurities  are  lo¬ 
cated  at  some  separation  from  this  layer  or  in  the 


bulk.  The  two-dimensional  electron  gas  can  be  made 
nondissipative  in  certain  ranges  of  large  applied  mag¬ 
netic  fields  at  low  temperatures,  when  these  conduc¬ 
tion  electrons  in  the  layer  are  in  the  integer  quantum 
Hall  effect  state.  Owing  to  this  property  and  also  larger 
qubit  separation  allowed,  we  consider  this  the  most 
promising  approach  and  focus  our  present  review  on 
such  systems. 

The  time  scale  of  the  qubit-qubit  interactions  will 
be  denoted  by  7jnt.  This  is  the  time  it  takes  to  accom¬ 
plish  a  two-qubit  quantum  gate,  such  as  CNOT  [2- 
5,57].  Typically  for  semiconductor  quantum  comput¬ 
ing  proposals,  Tint  <  Text,  and  in  fact  the  case  with 
Tint  <£  Text  has  some  advantages  because  one  can  use 
several  fast  single-spin  flips  to  effectively  switch  in¬ 
teractions  of  some  qubits  off  over  the  gate  cycle.  An¬ 
other  approach  to  controlling  (on/off)  of  the  two-qubit 
interactions  is  by  gates,  see  Fig.  2,  which  affect  the 
two-dimensional  electron  gas  and  the  localized  elec¬ 
tron  wavefiinctions. 

However,  the  same  conduction  electrons  that  pro¬ 
vide  the  qubit-qubit  interactions,  also  expose  the 
qubits  to  the  environment,  causing  relaxation  and  de¬ 
coherence.  Other  interactions  will  also  be  present,  that 
play  no  role  in  the  useful  quantum-computing  dynam¬ 
ics  but  contribute  to  these  undesirable  processes.  Re¬ 
laxation  and  decoherence,  and  their  associated  time 
scales,  are  addressed  in  the  next  section. 


3.  Time  scales  of  relaxation  and  decoherence 

The  processes  of  relaxation  and  decoherence  con¬ 
sidered  here  [21-28]  are  associated  with  the  dynamics 
of  a  small,  few-qubit  quantum  system  as  it  interacts 
with  the  environment.  Ultimately,  for  a  large,  multi¬ 
qubit  system,  many-body  quantum  chaos-like  behav¬ 
ior  must  also  be  accounted  for,  and  some  advances  in 
model  system  studies  have  been  reported  recently  [5, 
58].  Our  discussion  here  will  be  for  the  few-qubit  case 
mostly  because  it  allows  more  system-specific  investi¬ 
gations  for  actual  quantum-computing  proposals. 

Dynamical  processes  that  are  unwanted  in  quantum 
computing,  because  they  result  from  the  environmen¬ 
tal  influences  rather  than  from  the  controlled  radiation 
pulses  and  gate  potentials,  can  proceed  on  various  time 
scales.  In  fact,  it  is  not  guaranteed  that  processes  of 
various  types,  relaxation/thermalization  vs.  decoher- 
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ence/dephasing,  can  even  be  unambiguously  distinctly 
identified. 

At  low  temperatures,  it  is  generally  hoped  that  ther- 
malization,  which  requires  transfer  of  energy,  slows 
down.  If  the  fastest  such  processes  proceed  on  time 
scales  of  order  T\ ,  then  this  time  increases  at  low  tem¬ 
peratures  because  there  are  less  excitations  (phonons, 
electron  gas  modes,  etc.)  to  couple  the  small  quantum 
system  to  the  rest  of  the  solid-state  host  material. 

On  the  other  hand,  processes  that  do  not  require 
flow  of  energy  to  or  from  the  environment,  can  still 
effect  the  phase  of  the  quantum-superposition  ampli¬ 
tudes  and  cause  decoherence.  These  processes  can 
thus  proceed  faster,  on  the  time  scale  72.  While  these 
comments  seem  to  suggest  that  T2  ^  T\,  there  is  no 
obvious  reason  to  have  generally  T2  <3C  T\  at  low  tem¬ 
peratures. 

However,  if  the  spectrum  of  the  dominant  excita¬ 
tions  mediating  the  qubit  coupling  (both  to  each  other 
and  to  the  host  material)  has  a  gap,  then  we  expect  that 
all  the  relaxation  and  decoherence  processes  will  be 
suppressed.  Furthermore,  the  suppression  of  the  relax¬ 
ation  will  be  exponential,  with  the  Boltzmann  factor 
for  that  energy  gap.  Then,  T2  T\  will  be  satisfied 
but  also,  more  importantly,  the  actual  values  of  both 
time  scales  will  be  inordinately  large.  This  was  found, 
theoretically  and  experimentally,  to  be  the  case  for 
the  integer-quantum-Hall-state  two-dimensional  elec¬ 
tron  gas  as  mediator  of  the  localized-spin  (nuclear, 
electronic)  coupling  in  semiconductor  heterostructures 
[27,28,59-63]. 

It  is  important  to  emphasize  that  relaxation  and  de¬ 
coherence  are  really  many-body  properties  of  the  sys¬ 
tem  plus  environment.  Entanglement  with  the  environ¬ 
ment  owing  to  the  unwanted  couplings  results  in  the 
small  quantum  system  having  no  pure  wavefunction 
even  if  initially  it  was  prepared  in  a  pure  state.  Instead, 
it  can  be  described  by  a  statistical  mixture  represented 
by  a  density  matrix,  once  the  environment  is  traced 
over. 

This  reduced  density  matrix  of  the  system  is  ex¬ 
pected  to  evolve  to  the  thermal  one  at  large  times.  The 
approach  to  the  thermal  density  matrix,  which  is  di¬ 
agonal  in  the  system-energy  basis,  defines  the  time 
scale  T\.  If  the  temperature  is  low  enough,  then  there 
is  the  expectation,  see  [25,26]  and  references  therein, 
that  for  some  intermediate  time  scales,  of  order  72, 
the  density  matrix  becomes  nearly-diagonal  in  a  basis 


Initialize 


Control  and 
evolve  in  time 
Text  y  Tin,  «  T2,  T\ 


Control  qubits:  Text 

&  interactions:  Tin,  (  >  Tex, ) 


Avoid  relaxation:  T\ 

&  decoherence:  Tz  (<T{) 


Measure 

Fig.  3.  Evaluation  of  quantum  computing  models.  One  of  the 
criteria  for  feasibility  of  quantum  computing  in  a  given  physical 
system  is  the  possibility  of  initialization  of  the  qubits  in  the 
desired  superposition  state.  Another  important  design  consideration 
is  control  of  qubit  states  and  of  their  interactions.  In  order  to 
implement  quantum  computing  effectively,  the  time  scales  for 
realization  of  single  and  two-qubit  logic  gates,  Text  and  Tint* 
respectively,  should  be  several  orders  of  magnitude  smaller  than 
the  time  scales  of  relaxation  and  decoherence,  T\  and  Ti-  The 
relationships  between  these  time  scales  are  further  explained  in  the 
text.  Finally,  efficient  and  reliable  measurement  of  the  output  state 
of  the  qubits  is  required  for  reading  off  the  result  of  the  computation 
and  presently  represents  a  formidable  experimental  challenge. 


which  is  determined  not  by  the  systems  Hamiltonian 
(energy),  but  by  the  interaction  operator  with  the  en¬ 
vironment.  This  latter  process  corresponds  to  loss  of 
quantum  coherence. 

As  emphasized  in  Fig.  3,  evaluation  of  a  quantum¬ 
computing  proposal  requires,  among  other  things,  es¬ 
tablishing  the  relation  Tex t,  Tint  <3C  72,  T\.  Owing  to 
calculational  difficulties,  the  single-qubit  times  7^2 
will  usually  be  used,  though,  as  mentioned  earlier, 
some  study  of  the  multi-qubit  “quantum  chaos”  effects 
may  be  required.  For  spin-qubit  quantum  computing 
in  semiconductor  heterostructures,  the  relation  is  typ¬ 
ically  Text  Tint  ^  72  <$C  T\y  so  the  issue  is  usually 
how  small  is  the  quality  ratio  Q  =  Tint/ 72. 

The  required  value  of  Q ,  needed  for  fault-tolerant 
quantum  error  correction,  depends  on  the  physical 
model  of  error  sources  and  can  be  as  small  as  Q  = 
10_6-10-4,  see  [15,18-20],  or  as  large  as  Q  —  1/2, 
see  [64].  For  the  systems  of  interest  to  us  here,  spin 
qubits  in  semiconductor  structures,  the  value  of  Q  = 
10-5  is  a  reasonable  working  estimate.  Thus,  we  seek 
systems/conditions  with  Tmi/T2  ^  10-5. 


4.  Results  for  nuclear-spin  qubits 

In  this  section  we  outline  results  for  models  of 
quantum  computing  with  nuclear  spins  as  qubits,  and 
with  coupling  mediated  by  the  two-dimensional  elec- 
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tron  gas  in  the  integer  quantum  Hall  effect  state  [27, 
28,30].  In  strong  magnetic  fields,  the  spatial  states  of 
the  electrons  confined  in  the  two-dimensional  layer  in 
which  the  qubits  are  placed,  see  Fig.  2,  are  quantized 
by  the  field  to  resemble  free-space  Landau  levels.  The 
lattice  potential  and  the  impurities  actually  cause  for¬ 
mation  of  narrow  bands  instead  of  the  sharp  levels, 
separated  by  localized  states.  As  a  result,  for  ranges  of 
magnetic  field,  the  localized  states  fill  up  while  the  ex¬ 
tended  states  resemble  completely  filled  integer  num¬ 
ber  of  Landau  levels.  These  states  are  further  Zeeman 
split  owing  to  the  electron  spin.  At  low  temperatures, 
one  can  find  field  values  such  that  only  one  Zeeman 
sublevel  is  completely  filled  in  the  ground  state. 

The  electronic  state  in  such  systems,  that  show 
the  quantum  Hall  effect  [31]  in  conductivity,  are 
highly  correlated  and  nondissipative.  If  nuclear  spins 
are  used  as  qubits,  i.e.  atoms  with  nuclear  spin  1/2 
are  sparsely  positioned  in  the  zero-nuclear  spin  host, 
such  as  the  zero-nuclear-spin  isotope  28  of  Si,  which 
constitutes  92%  of  natural  silicone,  then  their  zero- 
temperature  relaxation  will  be  significantly  slowed 
down:  experimentally,  7j  ~  103  sec  [62]. 

Localized  spins,  both  nuclear  and  electronic,  inter¬ 
act  by  exchanges  of  spin  excitons — spin  waves  con¬ 
sisting  of  a  superposition  of  bound  electron-hole  pair 
states.  The  spectrum  of  these  excitations  [65,66],  ob¬ 
served  experimentally  in  [67],  has  a  gap  correspond¬ 
ing  to  the  Zeeman  splitting.  This  gap  is  the  cause  of 
slow  relaxation  and  decoherence.  The  exchange  of  vir¬ 
tual  spin  excitons  mediates  the  qubit-qubit  interaction 
and  also,  via  scattering  of  virtual  excitons  from  im¬ 
purity  potentials,  relaxation  and  decoherence  of  single 
qubits. 

The  original  proposal  to  use  nuclear  spin  qubits  di¬ 
rectly  coupled  by  the  two-dimensional  electron  gas 
[30],  required  positioning  the  qubits  at  distances  com¬ 
parable  to  several  magnetic  lengths.  The  latter  is  of  or¬ 
der  10  nm  for  magnetic  fields  of  several  Tesla.  The 
qubit-qubit  interaction  decays  exponentially  on  this 
length  scale.  Recently,  we  proposed  a  new  improved 
model  [28]  in  which  the  qubit  interactions  are  me¬ 
diated  via  coupling  of  the  two-dimensional  electron 
gas  to  the  outer  impurity  electrons.  This  applies  if  the 
atoms,  whose  nuclear  spins  are  the  qubits,  are  single- 
electron  donors  such  as  the  isotope  31  of  P.  These 
phosphorous  impurities  were  originally  utilized  in  the 
model  of  Kane  [32]  where  they  must  be  actually  posi¬ 


tioned  at  separations  of  about  4  nm  for  the  wavefunc- 
tions  of  the  outer  electrons,  which  are  bound  at  low 
temperatures,  to  overlap  significantly. 

In  our  new  improved  model  [28],  with  nuclear 
spins  coupling  to  the  outer  bound  electrons  which,  in 
turn,  interact  via  the  two-dimensional  electron  gas,  the 
interaction  turned  out  to  be  of  a  much  longer  range  as 
compared  to  the  model  of  [32]:  the  qubit  separation 
can  be  of  order  100  nm.  Another  advantage  is  that 
gate  control  of  the  individual  qubits  and  of  qubit- 
qubit  interactions  is  possible.  We  have  carried  out 
extensive  perturbative  many-body  calculations  [27,28, 
30,68]  allowing  estimation  of  and  Ti  for  both 
the  original  quantum-computing  proposal  [30]  and  its 
improved  version  [28],  where  the  main  improvement 
is  in  the  possibility  of  the  gate  control  along  the 
lines  of  [32].  The  “clock  speed”  of  the  improved 
model  is  also  faster  by  about  two  orders  of  magnitude. 
The  technical  details  of  these  rather  cumbersome 
calculations  are  available  in  the  literature  and  will  not 
be  reviewed  here. 

The  results  are  summarized  in  Table  1.  We  show 
estimates  of  all  four  relevant  time  scales  for  the  two 
models  introduced  earlier.  The  “original”  model  [30] 
corresponds  to  nuclear  spins  1  /2  introduced  at  qubits 
in  atoms  without  an  outer  loosely  bound  electron. 
The  “improved”  model  corresponds  to  the  case  when 
the  outer  electron  is  present  and  its  interaction  with 
the  nuclear  spin  and  the  two-dimensional  electron  gas 
dominates  the  dynamics. 

The  data  shown  in  Table  1  were  obtained  assuming 
typical  parameters  for  the  standard  heterojunctions 
utilized  in  quantum-Hall-effect  experiments  today, 
and  qubit  separation  of  65  nm.  Thus,  the  parameter 
values  taken  [28,30]  were  more  appropriate  for  the 
GaAs  system  than  for  Si,  even  though  the  main 
isotopes  of  gallium  and  arsenic  have  nuclear  spin 
3/2  and  cannot  serve  as  spin-zero  hosts.  The  reason 

Table  1 

Time  scales  of  the  qubit  dynamics  for  the  original  [30]  and  improved 
[28]  versions  of  the  nuclear  spin  quantum  computer  with  interac¬ 
tions  mediated  by  the  two-dimensional  electron  gas 


The  original  model 

The  improved  model 

Text 

0(1 0“5)  sec 

O(10-5)  sec 

^int 

0(1)  sec 

O(10-2)  sec 

h 

O(103)  sec 

0(10)  sec 

t2 

0(10)  sec 

OOO"1)  sec 
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for  using  these  values  has  been  that  experimental 
verification  of  some  of  the  numbers  might  be  possible 
in  the  available  materials  before  cleaner  and  different 
composition  materials  needed  for  quantum  computing 
are  produced. 

Our  estimates,  see  Table  1,  indicate  that  the  quality 
factor  Q  =  10-5  is  not  obtained  for  the  present 
system.  Actually,  no  quantum  computing  proposal 
to  date,  scalable  by  other  criteria,  satisfies  the  10-5 
quality-factor  criterion.  The  values  range  from  10-1  to 
1 0-2.  The  resolution  could  come  from  development  of 
better  error-correction  algorithms  or  from  improving 
the  physical  system  to  obtain  a  better  quality  factor. 
In  our  estimation  of  the  decoherence  time  scale, 
we  used  parameters  typical  of  a  standard,  “dirty” 
heterostructure  with  large  spatial  fluctuations  of  the 
impurity  potential.  These  heterostructures  have  been 
suitable  for  standard  experiments  because  they  provide 
wider  quantum-Hall  plateaus,  i.e.  ranges  of  magnetic 
field  for  which  all  the  extended  states  of  a  Zeeman 
sublevel  are  filled.  Much  cleaner,  ultra-high  mobility 
structures  can  be  obtained  by  placing  the  ionized 
impurity  layer  at  a  larger  distance  from  the  two- 
dimensional  gas  or  by  injecting  conduction  electrons 
into  the  heterostructure  by  other  means.  Thus,  our 
quantum-computing  proposals  [28,30]  are  unique  not 
only  in  the  large  qubit  separation  allowed  but  also 
in  that  there  is  a  clear  direction  of  exploration  to 
allow  physical,  rather  than  algorithmic,  resolution  of 
the  quality  factor  problem.  This  possibility  should  be 
further  explored  both  experimentally  and  theoretically. 
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Abstract 

Recent  theoretical  results  suggest  that  an  array  of  quantum  information  processors  communicating  via  classical  channels 
can  be  used  to  solve  fluid  dynamics  problems.  Quantum  lattice-gas  algorithms  (QLGA)  running  on  such  architectures  have 
been  shown  to  solve  the  diffusion  equation  and  the  nonlinear  Burgers  equations.  In  this  report,  we  describe  progress  towards 
an  ensemble  nuclear  magnetic  resonance  (NMR)  implementation  of  a  QLGA  that  solves  the  diffusion  equation.  The  methods 
rely  on  NMR  techniques  to  encode  an  initial  mass  density  into  an  ensemble  of  two-qubit  quantum  information  processors. 
Using  standard  pulse  techniques,  the  mass  density  can  then  manipulated  and  evolved  through  the  steps  of  the  algorithm.  We 
provide  the  experimental  results  of  our  first  attempt  to  realize  the  NMR  implementation.  The  results  qualitatively  follow  the 
ideal  simulation,  but  the  observed  implementation  errors  highlight  the  need  for  improved  control.  ©  2002  Elsevier  Science  B.V. 
All  rights  reserved. 


1.  Introduction 

The  field  of  quantum  information  processing  (QIP) 
has  made  steady  progress  over  the  past  decade,  driven 
in  part  by  the  realization  that  some  quantum  algo¬ 
rithms  offer  a  computational  advantage  over  the  best 
known  classical  counterparts  [1].  To  reach  a  practical 
improvement,  however,  most  quantum  algorithms  re¬ 
quire  a  large  number  of  qubits  coupled  quantum  me¬ 
chanically,  making  physical  implementation  difficult. 
Recently,  however,  it  has  been  suggested  that  some 
interesting  problems  might  be  solvable  by  a  hybrid 
classical-quantum  device  defined  as  a  type-II  quan¬ 
tum  computer  [2].  A  type-II  quantum  computer  is  es¬ 
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sentially  a  parallel  lattice  of  small  quantum  informa¬ 
tion  processors  that  share  information  through  clas¬ 
sical  channels.  Such  a  device  offers  the  experimental 
simplification  that  quantum  coherences  need  only  be 
maintained  locally  within  each  small  quantum  proces¬ 
sor.  Using  this  architecture,  it  might  be  possible  to 
increase  the  range  of  problems  that  small  quantum 
processors  can  tackle  by  classically  stringing  many  of 
them  together.  A  type-II  quantum  computer  may  thus 
serve  as  an  intermediate  architecture  between  few- 
qubit  and  large-scale  quantum  computers. 

In  this  report,  we  explore  the  experimental  aspects 
of  building  a  type-II  quantum  computer  using  nuclear 
magnetic  resonance  (NMR)  techniques.  Quantum  in¬ 
formation  processing  by  NMR  usually  employs  a  liq¬ 
uid  sample  of  molecules  containing  spin-^  nuclei  that 
is  subjected  to  a  strong  magnetic  field  [3].  A  typi- 
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cal  field  Bo  of  MO  T  creates  an  energy  difference 
A  E  between  the  aligned  and  anti-aligned  spin  states 
that  drives  the  system  to  an  equilibrium  state  with  net 
magnetization.  At  room  temperature,  A  E/kT  is  about 
10-5,  so  that  the  net  magnetization  is  relatively  small, 
but,  given  the  large  number  of  molecules  in  the  sam¬ 
ple  (~  1 0 1 8),  it  is  still  easily  detectable.  The  entire  spin 
ensemble  is  accurately  described  by  a  reduced  density 
matrix  of  only  the  intramolecular  spin  degrees  of  free¬ 
dom.  The  ensemble  nature  of  the  NMR  sample  thus 
makes  it  inherently  applicable  to  parallel  computation. 
A  type-11  architecture  can  be  mapped  onto  an  NMR 
sample  by  creating  a  correspondence  between  the  sites 
of  the  lattice  and  spatially  distinct  spin  ensembles.  Us¬ 
ing  magnetic  field  gradients  and  radiofrequency  (RF) 
pulses,  information  in  the  lattice  can  be  encoded,  ma¬ 
nipulated,  and  read  out.  As  a  first  test  of  the  NMR  im¬ 
plementation,  we  chose  a  basic  quantum  lattice  gas  al¬ 
gorithm  (QLGA)  that  solves  diffusive  dynamics  in  one 
dimension. 


2.  Lattice-gas  system 


The  diffusion  of  a  mass  density  p  is  governed  by 


dt  dt2  ' 


(i) 


The  above  equation  corresponds  to  the  macroscopic 
effective  field  theory  result.  Its  relation  to  the  lattice- 
gas  dynamics  is  seen  by  breaking  space  into  an  array 
of  lattice  sites  with  occupation  probabilities  assigned 
to  each  site  [4,5].  The  ensemble  average  mesoscopic 
dynamics  are  controlled  by  the  transport  equations  [2] 


f\(nym  -F  1)  =  f\(n,m)  +  \[f2(n,m)  -  f\(n,m)], 

(2) 

fl (h >  m  +  1)  =  /2(n,m)  -  ±[/2(n,  w)  -  /i(n,m)], 

(3) 


where  f\  and  /2  represent  occupation  probabilities 
and  the  bracketed  terms  represent  a  collision  operator. 
The  number  density  p  is  the  sum  of  f\  and  /2. 
The  indices  n  and  m  correspond  to  lattice  site  and 
time  step,  respectively.  The  connection  between  the 
diffusion  equation  and  the  transport  equations  may 
be  seen  by  taking  the  Chapman-Enskog  expansion 
of  the  lattice  Boltzmann  equation  written  in  terms  of 
occupation  probabilities. 


3.  Quantum  lattice-gas  algorithm 


The  quantum  lattice-gas  implementation  relies  on 
mapping  each  initial  occupation  probability  f\  and 
/2  into  the  corresponding  single-particle  states  of  two 
quantum  bits, 

|<7l,2(«,  m))  =  v//i,2(n,w)  |1)  +  v/1  -  f\,2(n,m)  |0), 

(4) 

where  \q\,i(n,m))  are  the  qubit  states  and  |0)  and  |1) 
correspond  to  the  eigenstates  of  a  two-level  system. 
The  resulting  two-qubit  wave  function  for  a  site 
becomes 


\'Hn,m))  =  y/fif2\U)  +  jM\  -/2)H0) 

+  Vo  -  /i)/2  |oi>  +  yd  -  /i)(i  -  7T)  loo)  (5) 


where  \r/r(n,m))  spans  the  Hilbert  space  of  two 
coupled  quantum  systems.  After  initialization,  the 
algorithm  calls  for  a  collision  operation 

\\lr'(n,m))  =  U\\//(n,m))  (6) 


that  is  carried  out  via  unitary  evolution  by  a  “square- 
root  of  swap”  gate  U .  The  gate  U  can  be  written  as 


U  = 


/I 

0 

0 

\o 


0  0\ 


1  -I-  ' 

2  +  2 

0 


2  +  2  0 

1  i 

2  "  2 

0 


0 

1/ 


(7) 


in  the  standard  basis.  The  next  step  in  the  computation 
requires  a  measurement  of  the  occupation  numbers 

=  (f\n,m)\n\'i\f\n,m)),  (8) 

where  the  number  operators  h\%i  are  defined  as 


p 

0 

0 

°\ 

p 

0 

0 

°\ 

0 

0 

0 

0 1 

1 0 

1 

0 

0 

n  i  = 

0 

0 
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0  • 

"2=  0 
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\0 

0 

0 

0/ 

\0 
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0/ 

(9) 


The  measured  occupation  numbers  f[  and  /2'  are 
then  streamed  to  the  nearest  lattice  sites  in  opposite 
directions,  as  given  by 


f\(n,m  -F  1)  =  f[(n  +  l,m),  (10) 

+  1)  =  7*2 (/i  -  1 ,  m).  (1 1) 

The  entire  diffusion  algorithm  can  be  summarized  in 
four  steps: 
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(1)  Initialization  of  occupation  probabilities  in  each 
spatially  distinct  site. 

(2)  Application  of  the  collision  operator,  U ,  at  all 
sites. 

(3)  Readout  of  the  expectation  values  of  the  number 
operators. 

(4)  Determination  of  the  new  occupation  probabilities 
by  streaming  to  nearest  neighbors. 


4.  NMR  implementation 

In  our  particular  test,  we  implemented  a  two-qubit 
diffusion  algorithm  using  a  solution  of  chloroform 
(13CHCl3)  where  the  hydrogen  and  the  labeled  carbon 
nuclei  served  as  qubits  1  and  2,  respectively.  Fig.  1 
shows  the  energy  level  diagram  of  the  spins  and  a 
picture  of  the  molecule.  As  shown  in  the  diagram,  the 
proton  splitting  is  four  times  larger  than  the  carbon 
splitting,  and  both  splittings  are  a  small  fraction  of  kT. 

4. 1.  Mapping  to  spin  ensembles 


The  first  step  in  creating  an  experiment  to  study 
the  implementation  of  type-11  quantum  computer  is 
to  define  a  mapping  of  the  theoretically  required 
quantum  states  to  a  real  physical  system.  In  the 
liquid-state  NMR  case,  the  required  quantum  states 
| \//(n,  m))  are  physically  encoded  onto  spin  ensembles 
described  by  density  matrices  o(n,m) 

\\lr(n,m))~*  a(n,m).  (12) 

However,  the  thermal  equilibrium  of  liquid-state  NMR 
systems  is  a  highly  mixed  state  that  is  not  immediately 
applicable  to  quantum  computing  experiments.  As  a 
result,  the  thermal  equilibrium  state  must  first  be  reset 
to  a  pseudopure  state  of  the  form  [6,7] 


Fig.  1.  The  picture  show  the  chloroform  molecule  and  the  nuclear 
spin  energy  level  diagram. 


<j(n ,  m)  =  1  —  e\ \J/(n,  m))(f(n,m)\.  (13) 

The  above  pseudopure  state  transforms  identically  to 
the  corresponding  pure  state  m)>.  Each  sub¬ 

ensemble  a(n,m)  is  in  turn  composed  of  a  large 
number  (MO18)  of  individual  molecules  distributed 
within  a  slice  of  a  cylindrical  sample.  More  formally, 
the  reduced  density  matrix  o(n,m)  at  a  site  is 


a(n,m)  =  Tr?^^r[z  -  Az(n  -  ^)] 

x  \4>{r,m))(<t>(r,m)\ 


(14) 


where  |0(r,/n)>  is  the  nuclear  spin  state  of  a  single 
molecule  located  at  position  r,  T (z)  is  the  “top  hat” 
function 


T(z)  = 


i.  izM, 

o.  \z\>{ 


(15) 


that  selects  the  relevant  spatial  slice  with  thickness  A z, 
and  Tr;  denotes  the  partial  trace  over  the  spatial  degree 
of  freedom.  The  variable  z  represents  the  correspond¬ 
ing  coordinate  of  the  vector  r.  The  states  a(n,  m)  re¬ 
quired  at  the  beginning  of  each  update  are  created 
by  applying  shaped  radiofrequency  (RF)  pulses  in  the 
presence  of  linear  magnetic  field  gradients.  This  step  is 
related  to  slice-selection  in  magnetic  resonance  imag¬ 
ing  (MRI).  Fig.  2  depicts  the  geometrical  arrangement 
of  the  slices  relative  to  the  gradient  and  RF  coils  in  the 
NMR  probe. 


Gradient  Coil 


RF  Coil 


Liquid  Sample  Separated 
into  Spatial  Nodes 


Fig.  2.  The  cylindrical  sample  of  chloroform  employed  in  this 
experiment  is  addressed  in  slices  by  the  combined  action  of 
magnetic  field  gradients  and  shaped  RF  pulses.  Each  slice  represents 
a  node  in  the  lattice  of  quantum  information  processors. 


342 


M.A.  Pravia  et  at.  /  Computer  Physics  Communications  146  (2002)  339-344 


Diffusion  Equation 


Fig.  3.  The  pulse  sequence  for  a  single  time  step  of  the  algorithm  begins  with  the  pseudopure  state  preparation.  Gradients  are  used  to  perform 
the  necessary  non-unitary  operations  required  for  equalizing  the  magnetization  of  the  two-spin  species  and  to  prepare  the  pseudo  pure  state.  The 
lattice  initialization  is  accomplished  by  applying  weak  RF  shapes  in  the  presence  of  a  magnetic  field  gradient  in  the  Z  direction.  A  decoupling 
sequence  prevents  the  scalar  coupling  from  interfering  with  the  initialization.  The  collision  operation  is  performed  by  a  sequence  of  coupling 
delay  and  strong  RF  pulses.  The  collision  pulse  sequence  is  applied  without  a  gradient  so  that  all  of  the  spins  in  the  lattice  feel  the  same 
operation.  Readouts  of  both  the  carbon  and  hydrogen  magnetizations  are  carried  out  on  the  hydrogen  channel  in  two  separate  experiments. 


4.2.  Control  and  measurement  of  spin  system 


In  the  absence  of  a  magnetic  field  gradient,  the 
Hamiltonian  of  the  spin  system  in  the  doubly-rotating 
frame  is 


H(t)  =  +  [u>'x(t)a}  +  1^(0^] 

+  [w2x(tt)ol  +  w2(t)cr2].  (16) 

The  first  term  denotes  the  scalar  interaction  between 
the  spins,  while  the  remaining  terms  are  the  externally- 
controlled  RF  Hamiltonian.  The  operators  of  the 
form  a x'%y%z  are  Pauli  spin  matrices  corresponding  to 
each  qubit,  and  the  scalar  coupling  Hamiltonian  is  a 
Kronecker  product  of  the  single-spin  operators.  The 
RF  part  of  the  Hamiltonian  generates  arbitrary  single 
spin  rotations  with  high  fidelity  when  the  nutation 
rates  w\'2y  are  much  stronger  than  7,  the  scalar 
coupling  constant. 

As  mentioned  before,  the  collision  operator  U  for 
the  diffusion  algorithm  is  the  square-root  of  swap  gate. 
The  unitary  operator  U  can  be  written  as 


U  =  exp 


(17) 


Hamiltonian  plus  appropriate  single-spin  rotations  [8]. 
The  operator  U  is  applied  to  all  the  sub-ensembles 
o(n,m)  such  that 

o' (n ,  m)  =  U a (n ,  m)U^ .  (18) 

The  final  steps  of  the  algorithm  are  to  read  the  oc¬ 
cupation  numbers  encoded  in  o'(n,m)  and  to  stream 
them  to  nearby  sites.  The  readout  is  accomplished  by 
noticing  that  Eq.  (8)  can  be  rewritten  in  terms  of  the 
z-Pauli  matrices  as 


f[2{n,m)  =  (y'(n, 


m) 


1  -  o 


1.2 


=  -{^\n.m)\al2\f(n,m))}  (19) 

using  the  fact  that  h\2  =  j(l  —  a}'2).  The  last 
equation  can  be  written  in  the  final  form 


(20) 


where  the  trace  has  been  replaced  by  the  z-magnetiza- 
tion  A/i  2.  The  z-magnetization  is  measured  in  NMR 
by  applying  a  “read”  7t/2  pulse  and  observing  the 
transverse  magnetization.  The  measured  values  /,'  2 
can  be  streamed  on  a  classical  computer  and  then 
reinitialized  onto  the  lattice. 


if  an  irrelevant  phase  is  ignored.  Written  in  this 
form,  it  is  clear  that  U  can  be  decomposed  into  the 
product  of  three  commuting  terms.  Each  term  can  be 
implemented  by  making  use  of  the  scalar  coupling 


4.3.  Pulse  sequence 

The  diagram  in  Fig.  3  shows  the  main  parts  of  a  sin¬ 
gle  time  step  of  the  NMR  implementation:  pseudopure 
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Simulation  of  Readouts 


Fig.  4.  The  plots  show  a  simulated  run  of  the  ideal  quantum  lattice  gas  algorithm  for  diffusion.  The  top  left  plot  contains  the  first  step,  followed 
to  the  right  and  then  down  the  rows  by  subsequent  time  steps.  The  dashed  lines  represent  the  initialized  occupation  numbers  f\  2  for  each  spin, 
while  the  solid  lines  represent  the  occupation  numbers  /[  2  present  after  the  collision.  The  jc-axis  labels  the  node  number  and  runs  from  1  to  16. 


Readouts  at  Each  Step 


Fig.  5.  The  experimental  results  for  the  corresponding  time  steps  of  the  simulations  from  Fig.  4.  Although  the  general  features  follow  the 
simulation,  the  experimental  results  are  not  of  high  fidelity  and  suggest  a  need  for  more  precise  control.  The  jc-axis  labels  the  observed  spectral 
frequency.  The  actual  nodes  used  in  the  experiment  reside  in  the  region  between  —200  and  200  Hz.  The  outlying  region  is  included  for  reference. 
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state  preparation,  lattice  initialization,  collision,  and 
readout.  The  top  two  lines  correspond  to  operations 
on  the  two  qubits  (H  and  C),  while  the  third  line  shows 
the  required  gradient  pulses.  The  pseudopure  state  was 
prepared  by  first  equalizing  the  magnetizations  of  the 
two  spins,  followed  by  a  pseudopure  state  creation  se¬ 
quence  [9].  The  starting  occupation  numbers  for  each 
time  step  were  then  encoded  using  weak  shaped  RF 
pulses  on  the  two  spins.  Because  the  RF  power  utilized 
was  weak  relative  to  the  gradient  strength,  the  shape  of 
the  pulse  was  determined  by  taking  the  Fourier  trans¬ 
form  of  the  desired  magnetization  profile  [10].  A  de¬ 
coupling  sequence  was  applied  simultaneously  with 
the  RF  shape,  to  average  out  the  effects  of  the  scalar 
coupling  on  the  RF  excitation. 

The  collision  operator  was  implemented  by  decom¬ 
posing  the  total  unitary  operator  into  sequences  of 
scalar  coupling  delays  and  RF  pulses.  The  readout 
was  performed  by  recording  the  spectra  in  the  pres¬ 
ence  of  a  weak  gradient.  The  classical  communica¬ 
tion  part  of  the  algorithm  was  absorbed  into  the  en¬ 
coding  operation  of  the  next  time  step.  A  linear  phase 
ramp  was  added  to  the  RF  shape,  effectively  shift¬ 
ing  the  frequency  of  the  excitation.  Since  the  data  on 
the  two  spins  was  to  be  shifted  in  different  directions, 
the  phase  ramps  for  the  two  RF  pulses  had  opposite 
slopes. 


5.  Results 

The  two  Figs.  4  and  5  show  the  results  of  pre¬ 
liminary  experiments  and,  for  comparison,  the  cor¬ 
responding  ideal  simulation  of  the  NMR  implemen¬ 
tation.  The  experiments  where  performed  using  16 
nodes  iterated  through  12  time  steps  of  the  algorithm. 
As  can  be  seen,  the  broad  features  of  the  diffusion 
can  be  seen,  but  large  errors  are  present  in  the  imple¬ 
mentation.  The  errors  are  caused  by  problems  in  the 
decoupling  sequence,  errors  in  the  Fourier  transform 
approximation,  and  other  experimental  imperfections. 
We  are  continuing  to  refine  the  experiments  and  we 
expect  to  correct  these  errors  in  the  near  future  with 
improved  results  to  be  published  in  a  subsequent  pa¬ 
per  [11]. 


6.  Conclusion 

Ensemble  NMR  techniques  have  been  success¬ 
fully  used  to  study  the  experimental  details  involved 
in  quantum  information  processing.  The  astronomi¬ 
cal  number  of  individual  quantum  systems  (~1018) 
present  in  typical  liquid-state  spin  ensembles  greatly 
facilitates  the  problem  of  measuring  spin  quantum  co¬ 
herences.  In  addition,  the  ensemble  nature  of  the  sys¬ 
tem  has  been  successfully  utilized  to  create  the  nec¬ 
essary  pseudo-pure  states  and  to  systematically  gen¬ 
erate  non-unitary  operations  over  the  ensemble  [12]. 
In  this  implementation,  we  again  exploit  the  ensemble 
nature,  but  this  time  as  a  means  of  realizing  a  lattice 
of  quantum  information  processors.  The  implementa¬ 
tion  combines  the  advantages  of  quantum  computation 
at  each  node  with  parallel  computation  throughout  the 
lattice.  The  large  size  of  the  NMR  ensemble  provides, 
in  principle,  sufficient  room  to  explore  large  lattices. 
Although  achieved  experimental  results  point  to  the 
need  for  better  control,  the  experiments  are  a  first  step 
towards  realizing  the  quantum  lattice  gas  algorithm  on 
a  NMR  quantum  information  processor. 
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